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ABSTRACT 

Wc report on a very large set of simulations of collisions between two main sequence 
(MS) stars. These computations were done with the "Smoothed Particle Hydrody- 
namics" method. Realistic stellar structure models for evolved MS stars were used. In 
order to sample an extended domain of initial parameters space (masses of the stars, 
relative velocity and impact parameter), more than 14000 simulations were carried 
out. We considered stellar masses ranging between 0.1 and 75 Mq and relative veloc- 
ities up to a few thousands kms~^. To hmit the computational burden, a resolution 
of 1000-32000 particles per star was used. The primary goal of this study was to 
build a complete database from which the result of any collision can be interpolated. 
This allows us to incorporate the effects of stellar collisions with an unprecedented 
level of realism into dynamical simulations of galactic nuclei and other dense stellar 
clusters. We make the data describing the initial condition and outcome (mass and 
energy loss, angle of deflection) of all our simulations available on the Internet. We find 
that the outcome of collisions depends sensitively on the stellar structure and that, 
in most cases, using polytropic models is inappropriate. Published fitting formulae for 
the collision outcomes, established from a limited set of collisions, prove of limited use 
because they do not allow robust extrapolation to other stellar structures or relative 
velocities. 

Key words: hydrodynamics - methods: numerical ~ stars: interior - galaxies: nuclei, 
star clusters - Galaxy: centre 



1 INTRODUCTION 

1.1 Stellar collisions in galactic nuclei 

In the recent years, the study of stellar collisions has received 
renewed interest from researchers studying the dynamics of 
dense stellar systems, either open/globular clusters_or galac- 
tic central regions (see the contributions in lSharal2002h . Our 
own motivation is to perform simulations of the long-term 
evolution of dense stellar systems, particularly galactic nu- 
clei, with a new Monte Carlo stellar dyna mics code which 
incorporate collisions as "micro-physics" iFreitag fc Benj 
12001 hll2002>) . 

Before going into a brief description of the astrophysi- 
cal motivations of these works, some clarification about the 
notion of "stellar collision" is called for. In this paper we 
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shall use this term to refer to a process during which two 
stars, previously unbound to each other, get so close that not 
only gravitational forces but also hydrodynamical ones come 
into play to determine the outcome of the interaction. So, 
strictly speaking, collisions are not only contact encounters 
but also comprise tidal interactions. However, for reasons to 
be exposed in Sec. 12.81 we restrict here to events leading to 
physical contact at first periastron passage. 

To assess the importance of collisions in a given astro- 
physical context, the key quantity to monitor is the collision 
time, which we define here as the average time needed for 
"test-star" 1 to experience a collision with any "field-star" 
2, 

Tcoii = {n2Svrciy^ (1) 
with S = TT^oii 1 H ^- 2 '- , (2) 

where 712 is the number density of the field stars, Wroi the 
relative velocity and dcoii the pericentre distance leading to 
a collision (dcoii = Ri + R2 for contact at first passage, ne- 
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Figure 1. Relaxation and collision times at 0.1 pc from a massive black hole in the centre of a galactic nucleus (inspired from a similar 
diagram bv lArabadii3ll997l) . We plot curves of iso-Ti-dax (dashes, blue in colour version) and iso-Tcoii (solid lines, magenta in colour 
version) in a plane parameterised by the stellar density at 0.1 pc and the mass of the central black hole. The right ordinate scale indicates 
the stellar velocity dispersion (Keplerian value, ay = \/ GA/bh /''')■ AH stars are assumed to be of solar type. The left shaded sector 
(cyan in colour version) corresponds to conditions for which both T^eiax ^■nd 'TtLoW ^•rfi larger than Hubble time, so that such nuclei are 
not expected to show signs of secular evolution. In the shaded lower right corner, the black hole does not dominate the kinematics and 
the effect of the cluster's self gravitation should be included in the computation of the characteristic time-scales. 

Large black dots show the estimated conditions for observed galactic nuclei. In most cases, the estimation of the stellar density at 0.1 pc 
requires important extrapolation from the data, as such a small radius is resolved only for a few galaxies of the local group (the Milky 
Way, M31 and M32). For th is extrapolation , we us e a power-law cusp of the form p(r) = Po(r'o /r)^. The va l ues of poi ^0 a^nd 7 are taken 
from lGebhardt et al.H2003h o r lFaber et al.l <1997ft . The densities for M31 and M 32 are fr o mlLauer et al.Nl998l) and the Milky Way's 
value frorr nc'enze^talTiMg^ . The values of Mbh are estimates by Magorr ian et al van der Marell Jl999ft or better constrained 

values gathered by Kormendy l2004t see http : //chandra. as .utexas . edu/'kormendy/bhsearch .html for these data and a list of original 
references). In some cases, po and ro have already been extrapolated from larger radii! Cases w ith an horizontal l ine connected to a 
second smaller dot are nuclei for which the slope 7 is observational compatible with 0, according to lFaber et al.l il997t) . The second point 
indicates the density value at 0.1 pc if constant p is assumed up to ro. 



Comprehensive set of main sequence stellar collisions 3 



giecting tidal deformation). 5* is the coUisional cross-section. 
At low relative velocity, it is greatly enhanced over the geo- 
metric value by gravitational attraction. This effect, dubbed 
"focusing" is expressed by the second term in the brack- 
ets of Eq. 121 In most astronomical contexts, the velocity 
dispersion is much smaller than the stellar escape veloc- 
ity K = ^/2GM7[RI (~ 620kms~^ for sun-like stars) and 
gravitational focusing dominates. In these cases, integrating 
over a Maxwellian distribu tion for relative velocities yields 
jBinnev fc Tremainelll987l Eq. 8-125): 

Tcoii ~ 7 X 10" yrs x (3) 

\Mq) Uc-V Vkms-i; 

In systems with Tcoii smaller than typical stellar ages, 
collisions have expectedly imprinted not only the stellar pop- 
ulation but also the global dynamical structure. Very high 
densities are necessary for such situations to take place but 
even when collisions occur at frequencies too low to be of 
dynamical relevance, they still can be of great astrophysi- 
cal interest per se because they are suspected to lead to the 
formation of unusual individ ual stellar objects, such as blue 
stragglers or stripped giants jDaviejll99a ISharalll99gl and 
references therein). Collisions are unimportant in the bulk 
of a galaxy; the probability for the sun to suffer a collision 
during its 10 Gyr main sequence life, amounts to no more 
than 10~^! Only in stellar clusters and galactic nuclei, is 
there a non- vanishing probability for at least some stars to 
experience collisions. For reviews about the role of collisions 
in var i ous en vironments, we refer to the various papers in 
|ShOTi(|200i). 

Among known stellar environments, galactic nuclei are 
those in which the most extreme values of the stellar density 
and velocity dispersion are attained. The best known case 
is our own Galaxy. Inside a sphere of radius 0.4 pc at the 
Galactic centre, the stellar density exceeds 4 x 10® Mqpc"'^ 
and a velocity dispersion of order 500kms~^ has been re- 
ported at a distance of 0.0 1 pc of the 2 — 3 x 10® M0 central 
black hole jGenzel et aljIlQQG. .20001) . Most other galactic 
nuclei are not resolved yet so we can only produce very un- 
certain estimates of Tcoii for these systems. Some of them are 
indicated in Fig. Q As a bias toward our own interests, we 
treat only the situation of a massive black hole dominating 
the kinematics of the surrounding stars. 

From this diagram, we see that there are very few galax- 
ies for which we can be certain that collisions played a im- 
portant dynamical role. Using Nuker model fits to represent 
the density profiles and the empirical relation between the 
mass of the central obj ect and the velocity d ispersion in the 
spheroidal com pon ent jTremaine et alJl2002l) as a proxy for 
the BH's mass. fYul i2003l) estimated the collision times for a 
series of observed galactic nuclei. She found only a few cases 
with Tcoii possibly shorter than the Hubble time and that, 
in present-day nuclei, collisions do not produce observable 
colour gradients in the stellar populations. It may be that 
the importance of these processes has been somewhat over- 
estimated in the past Ivan den Bcrgh,.1965; .Sanders 1970^1) . 

The centre of the Gal axy is a particularly complex and 
fascinating environment ijOenzel et alJ l2003t lOhez et alJ 
1200,4 ISchodel et all l2003^ . The "SO" stars orbiting the 
3 - 4 X 10® Mq BH Sgr A* at distances smaller than 0.04 pc 



seem to^^oi^ft^^^^iti^nasse^^^fias^^^^^^Ghez 
et al.TgOO^^Recen^teUa^form^oi^^tifi im- 
possible and scenarios to bring them from a few pc away in 
less than their sh ort lifetime require considerable fine tuning 
iKim et alj|2004l and references therein). Consequently, it is 
tempting to hypothesise they were created in a sequen ce of 
mergers of older, hghter MS stars dCenzel et al.ll2003l) . Us- 
ing simple Fokkcr-Planck modelling (not including a central 
BH'l. iLec 1,1994 199fift concluded that mergers can not ac- 
count for the formation of the massive stars found near the 
centre. On the other hand, whether collisions are responsi- 
ble for the observed relative depletio n of red giants at the 
Gala ctic centre is still a d ebated issue ( Gerhar dll994l: Davies 
et al. ll998HAlexandeill999.: ,Bailcv fc Davie jll999^ . Clearly, 
more detailed stellar dynamical models, that take into ac- 
count the presence of the central BH and include a realistic 
treatment of collisions and stellar evolution are called for to 
establish the role of collisions in the MW central cluster. 

There are however strong theoretical motivations to be- 
lieve that stellar encounters may have taken place in large 
numbers in the past evolution in many galactic centres with 
sufficiently high stellar densities. The main reason is the 
presence of massive compact dark objects in the centre 
of many, if not most, galaxies. These mass concentrations 
are most probably supe rmassive black holes (SMBHs) with 
mass e s lO'^-Sx 10^ Mr:. jKormendv fc Richstonell995HBarthl 
l2004HBarth. Ho. Rutledeej fc Sargend 120041: Fe rrarese & 
Ford l2004HKormendvll2004t iPinknev et aljl2003ft . From a 
series of works publishe d in the 70 's llPeeblejll97l Shapiro 
fc Light man 1976: BahcaU fc Woljll976lll977l: Dokuchaev & 
Ozernoi' IT977albt lYoung|ll977bl among others), it is known 
that a SMBH-surrounding stellar system whose long-term 
evolution is driven by 2-body gravitational encounters will 
relax to a density profile, close to a power-law p(r) oc r 
which yields a constant ffux of stars toward the centre where 
they are destroyed either by tidal disruptions or energetic 
collisions. If all stars have the same mass, the exponent is 
7 = 1.75. In the innermost regions of such a cusp, a high 
collision rate is expected. But the collisions themselves could 
act as a feed-back mechanism on the evolution of the stellar 
system and the growth of the black hole so that the actual 
formation of relaxa tional cusp is questionable. From analyt- 
ical considerations. iFrankI ^1978^ concludes that collisions in 
the cusp are never of importance, when compared to tidal 
disruptions, but this statement is seriously challenged by 
other studies and, in particular, more recent numerical sim- 



ulations lYoung et al.ll97'^ 


Young|l977aHDuncan fc Shanird 


ll983HDavid et alJll987aIbl: 


Murnhv. Cohn. fc Durisenlll99ll: 


lRauchlll999l). Unfortunatelv. the discussion of the contri- 



bution of various dynamical processes to the evolution of 
galactic nuclei has been blurred by uncertainties about the 
precise outcome of these individual processes. For instance 
the amount of gas that is accret ed by the SMBH following a 
tidal disruption is still debated llAval et al.ll2000l and refer- 
ences therein). As for stellar collisions, most previous works 
relied on quite unrealistic prescriptions, like complete de- 
struction fYoung et al.lll977t lYounelll977at iMcMillan et alJ 
[1981: Duncan fc ShapirJl983ft o r on a semi-analytical recipe 
proposed bv ISpi^e^^Sash^ il966l hereafter SS66) that 
complet ely neglects the real hydrodynamica l nature of the 
process (ISanderslll970aHDayid et aljlll987al5 iMurphv et alJ 
I1991I) . The work of iRauchI Jl999ft is a noticeable exception; 
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he used the results of a set of hydrodynamical simulations 
of stellar collisions by Davies to derive fitting formulae for 
the quantitative outcome of these events. The present work 
originated in our wish to get rid of these annoying uncer- 
tainties about the role of collisions in dynamica l simulations 
of galactic nuclei (iFreitag fc Benj|2001all2002l) . 

Many of the papers we have just cited were not only 
concerned with the past evolution of galactic nuclei but 
also (or mainly) with scenarios to feed SMBH and pro- 
vide quasars' luminosities. Gas-dynamical processes are now 
favoured candidates for the fuelling of active galactic nuclei 
(AGN) and th e dense cluster hypothesis seems so mewhat 
out-of-fashion JShlosman et al.lll99(]l: ICombedlioOll and ref- 
erences therein). On the other hand, AGN models have been 
proposed in which large luminosities in electro-magnetic ra- 
diation and/or relativistic particles are emitted by the hot 
gas clouds created by very e nergetic stellar collisions. First 
propositions along that line JWoltieJIlQel ISanderj|l970d) 
postulated that the stars' velociti es were due to the cluster's 
self gravitv. More recent models llKeena Dokuchaev 
et al1l993HCourvoisier et alllQQfitlTorricelli-Ciamponi et alj 
l200(j) invoke a SMBH to provide velocity dispersions rang- 
ing from a few lO'^kms"^ to a few lO'^kms^^. These non- 
standard AGN models may be successful in reproducing ob- 
served luminosity- variability relations that are otherwise dif- 
ficult to explain, but they should be re-examined in the light 
of a more refined treatment of stellar collisions and stellar 
dynamics. A third possibility for collisions to contribute di- 
rectly to the luminosity of AGN is to boost the rate of super- 
novae through creation of mass ive stars by mergers (IColgatd 
ll967tlShields fc WheeleJll978h . 

Finally, even though they are not the dominating lumi- 
nosity source in AGNs, stellar collisions may be responsible 
for the formation of massive black holes in dense galactic 
nuclei, eithB^bvrm^wa^ne^in^|Sande^^^70a^(3uinlan 



Freitag, & Giirkan |2004 


Giirkan. Freitag. & Rasia .2o61 


iFreitag. Giirkan. & Rasid 


2004blal |2005a|) or bv creating a 



massi ve gas cloud that subsequ ently evolves to a black hole 
iSpitzcr fc Saslaw..l96&: .Bcgclrnan fc Reejll978t Langbein 
et al.ll990l). 



1.2 Previous simulations of collisions between 
main sequence stars. 

Table lists the previous computations of high-velocity 
collisions between MS stars. We only mention "realistic", 
multi-dimensional hydrodynamical simulations. This ex- 
cludes early calculations that were based either on semi- 



analytical methods or on one-dimensional numeri- 

cal schemes (lMathij|l967t DeYoung|l968^ . Such approaches 



were clearly over-simplifications in which the real 3D hy- 
drodynamical nature of the problem was not properly ac- 
counted for. The importance of these "pre-hydrodynamics" 
works should not be underestimated, however. For instance, 
even though it was always deemed too simplistic to yield 
better than order-of-magnitude estimates, the lSS6(3 method 
had been adapted and used in a few key simulation works. 
We postpone a presentation of this "historical" method to 
Sec. I3.2l where we compare our results to predictions of this 
approach. 

With the historical exception of ISeidl fc CameronI 



1(1972), s-l' cited works were realised using SPH (Smoothed 
Particle Hydrodynamics). When featured with a tree to 
compute gravitation, SPH is a grid-less method which can 
cope with any asymmetrical three dimensional geometry. It 
ignores void spaces completely, imposes no physical limits 
beyond which matter cannot be tracked, does not come into 
trouble with large dynamic range as long as variable smooth- 
ing lengths are implemented. SPH is better suited to highly 
dynamical problems than to near-equilibrium configurations 
JSteinmetz fc Mulleiliggl) . For all these reasons, SPH is par- 
ticularly well suited to the simulation of stellar collisions. 

From Table Q it is clear that the study of the outcome 
of high- velocity collisions has not attracted much interest in 
the last few years, in contrast to parabolic encounters (Lom- 
bardi et al. Il996t ISills fc Lombardilll997t ISills et alj|200ll 
I2OO2I among others). As a consequence, the resolutions used 
seem very modest , by present-day standards; for instance, 
ISills et al.l l|2003) present a parabolic collision simulated 
with 10^ particles. Obviously, the simulations presented in 
this work do not correspond to a break-through in terms 
of resolution. This reflects the fact that most computations 
were realised a few years ago, when computing power was 
more limited and, most importantly, that we had to cover 
a huge parameter-space, requiring more than 10 000 simula- 
tions (see Sec. 12.71 . This sheer quantity, combined with the 
use of realistic stellar models instead of polytropes, represent 
the main improvements over previous works. 

For simpl icity, in the remaini ng of this ar ticle, we refer 
to the work of lBenz fc ffilS Jl987l) aslBHSTl tolBenz fc Hilh 
il992l) as BH92. tolLai Rasio. fc Shapirol lll993D as lLRS9£ 
and to lRauchI lll999h as lR99l For a more comprehensive list 
of references on simulations of all kind of stellar collisions, 
see the web site maintained by MF in the framework of the 
"MODEST" collaboration^ . 



1.3 Collisions with non-main-sequence stars 

In this work, we only treat collisions between two main- 
sequence (MS) stars. The motivations for this choice was to 
keep the number and variety of collisions to consider at a 
manageable level and that the present version of our Monte 
Carlo code only includes simplified stellar evolution which 
skips over the giant phase and turns MS stars directly into 
remnants. However, in a real stellar system, MS-MS encoun- 
ters may not dominate the global collision rate. Indeed, a 
given star of mass > 1 Mq may have a smaller probability 
for colliding with another star during its MS life than dur- 
ing its red giant (RG) phase despite t he latter being about 
10 times shorter (iBressan et al.lll993l for instance). This is 
made very clear by integrating the coUisional cross-section 
over the lifetime of the star, as we did in Fig. |21 In many 
cases, the collision probability during the red giant phase 
exceeds its MS counterpart for high relative velocities. RG- 
RG collisions are less likely than RG-MS events. Indeed, the 
ratio of probabilities can be estimated as follows 



ftlG-MS "-MS 

-Prg-rg n-fiG 



Rrg 



(2i?B 



: 0.25 



Tms 
2rg 



10. 



1 "MODEST" stands for Modelling DEnsc STcUar systems, see 
http://www.maiiybody.org/modest/ For the collision "working 
group", go to http://www.manybody.org/modest/WG/wg4.htmlj 
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Table 1. Hydrodynamical simulations of high- velocity collisions between MS stars in the literature. 



Reference 


Abbrev. 


Stellar models 




q = M-1/M2 




M—R relation 


Method 


Seidl & Cameron ('1972) 




polytropes n = 


3 


1 


0, 1.6, 3.2 




Head-on, 2D finite diff. 


Benz & Hills ("1987) 


BH87 


polytropes n = 


1.5 


1 


0-2.33 




SPH 1000 part. 


Benz & Hills (1992) 


BH92 


polytropes n = 


1.5 


0.2 


0-1. 5('^) 


R, oc Mj?-**^ 


SPH 7000 part. 


Lai et al. (199^^ 


LRS93 


polytropes n = 


3(d) 


0.1-1 


0.2-3.8 


R* oc Mj?-** 


SPH 8000 part. 


Davies (Ranch 1999") 


R99 


polytropes n = 


3 


0.25, 0.5, 1 


1-25 


i?, oc A/, 


SPH ~ 40 000 part. 


This work 




realistic 




0.0013-1 


0.1-30 


realistic 


SPH 4000-36 000 part. 



See symbols definition in Sec. 12.11 
''■'Up to 5 for head-on collisions. 

Fitting formulae are given. 
(■^^Eddington models. 
(°) Results only given as fitting formulae. 




Figure 2. Cumulative collision probability (normalised to 1) integrated over the lifetime of three stellar models. The second star is 
assumed to be sun-like. In each case, three velocity regimes are considered. At low relative velocities (V^^ <C V*), the collisional cross- 
section scales like i?* + Rq (solid lines); at very high velocities (V^^ ^ Vt), we recover the geometrical cross section, oc + R©)^ 
(dashed lines) . We also plot the case for a relative velocity of 200 km s~ ^ , an intermediate value typical for galactic nuclei. The end of the 
M S phase corresponds to the point where the slope of the curves increases suddenly. Stellar evolution models are from the compilation 
of lLeieune &: Schaererl J200ll) . available on-line at |http : //vizier ■ cf a .harvard. edu/vlz-bin/VlzieR?-source=VI/102| 



Although probably more common than MS-MS encounters, 
RG-MS collisions may not be more important. RG envelopes 
have very low densities so only little mass is lost in most 
cases and the RG recovers its appearance. At relative veloc- 
ities found in galactic nuclei, the MS star cannot be captured 
unless it is aimedjM^l^^iregU^^^tii^^^^^en^^^Bailey 
& Davies 1999). Furthermore, as giants will loose their en- 
velope anyway through winds and a planetary nebula or SN 
phase, collision with giants will probably make little differ- 
ence as far as the feeding of a central SMBH is concerned. 

Due to mass segregation in clusters and nuclei, collisions 
between compact remnants (CRs) and MS (or RG) stars are 
probably much more common than the small CR fractional 
number would suggest. For instance, the innermost 0.1 pc 
of the Sgr A* cluster is likely dominated by invisible stellar 
BHs (Miralda-Escudc fc Gould.200(]i) which may collisionally 
destroy MS and RG stars. CR-MS and CR-RG coUisions 
may also be of great interest as a way to produce exotic 
objects, such as cataclysmic variables. Unfortunately, due 
to the high dynamical range involved, the hydrodynamical 



simulation of these events is challenging and comprehensive 
predictions for their outcome are still lacking. 

Now that the astrophysical stage is set, we can proceed 
with a description of our simulation work. In Sec. |21 we de- 
scribe the choice and setting of initial conditions and present 
the numerical methods we use to compute and analyse stel- 
lar collisions. In Sec 13 results are reported and we explain 
how to exploit them in stellar dynamical simulations. Fi- 
nally, in Sec|ll some general conclusions and a discussion of 
further work to be done are presented. 



2 DESCRIPTION OF THE APPROACH 

2.1 Definitions, basic formulae and units 

As some quantities will be referred to again and again, we 
find it useful to define them once for all at the beginning 
of this article. Collisions between two main sequence (MS) 
stars are considered. In the centre of mass frame, the colli- 
sion is completely determined by 4 quantities: the masses Mi 
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Figure 3. Sketch of a gravitational 2-body hyperbolic encounter 
in the centre of mass reference frame. 



and M2 (in our work, we made the unconventional notation 
choice: Mi ^ M2), the impact parameter b (see Fig. O and 
the relative velocity at infinite separation, V^'^. The stellar 
radii are Ri and R2. Instead of applying a simple but un- 
realistic power-law mass-radius relation, the values for the 
radii are taken from the stellar models described in Sec. 12.21 
We shall often refer to the situation of a 2 point-mass 
hyperbolic encounter where all finite-size (hydrodynamical) 
effects are neglected. In this case, we define the periastron 
distance, 



{Ri + R2) 



1 + \/l +462^)4 



(4) 



with b = 6/(i?i + R2) and v = K^/K (see Eq.lSJ. When 
gravitational focusing is important, dmin is a more conve- 
nient parameter than b? Ignoring tidal effects such as de- 
formations and trajectory modification until strong hydro- 
dynamical interactions begin, we assume that only collisions 
with dmin < Ri + R2 lead to contact between stars at the 
first periastron passage. 

In case both stars survive the encounter and are left un- 
bound to each other, we define a coUisional defiection angle 
Scoii- This is the angle between the direction of the initial 
relative velocity (at infinite separation) with the direction of 
the final relative velocity (at 00). To assess the importance 
of finite-size hydrodynamical effects, it is useful to compare 
dcou with the value for a Keplerian hyperbolic orbit Sgrav, 



tan 



G(A/i + A/2) 



(5) 



A natural velocity scale for collisions is the relative ve- 
locity at contact for stars initially at rest at infinity. 



K = 



2G(Afi + M2) 
R1+R2 



(6) 



The structure of MS stars with A/* > 0.5 Mq is very 
concentrated (see Fig 0] and the appendix) and the total 
radius is not a good indicator of the extension of the stellar 
matter. It is thus often useful to normalise quantities with 
reference to the half- mass radius i??', i.e. the radius of a 
sphere that contains half the stellar mass. These radii can 
be read from the 50 % curve in Fig.gf. We can then define 




Figure 4. Mass-Radius relation for the stellar models. The thick 
line shows the total radius as a function of the mass of the star. 
Thin lines are for inner Lagrangian radii containing 25 to 95 % 
of the mass. Data for M<. ^ 0.4 M(7) (dots) is from realistic stel- 
lar structure models jSchaller et al.lll992t iMevnet et al.lll994 
ICharbonnel et al.lll99gh . Below 0.4 Mf;^. a polytropic (n = 1.5) 
structure is applied with a power-law M-R relation extrapolated 
from higher masses. The run of various Lagrangian radii makes it 
obvious that stellar models for different stellar masses cannot be 
deduced from each other through any homologous scaling. The 
"real" mass, determined through relation |2l is given in abscissa, 
not the ZAMS value. 



a "half-mass velocity" scale through 

(h)^ l 2G{Mi+M2) 

" ' R(h) , „(h) ■ 



(7) 



This quantity gives a better idea of the relative velocity when 
strong hydrodynamical effects begin to play an important 
role. Note that we use total masses in this definition. We 
often normalise initial parameters by these half-mass quan- 
tities so handy definitions are 



A : 



and u 



(h) ■ 



(8) 



Rq 

Mo 
V© 



Typical scales are set by "solar units", i.e.: 

7 X 10** m, 
2 X 10^" kg. 



:= (GMq/R©)'/^ = 436.5 kms"\ 
To := (RoV(GMo))'^' = 1604 s. 
These values are also referred to as "code units" 



^ Furthermore, b is not defined for parabolic encounter whereas 
djnin still is. 

^ We could also use radii at 75, 90 or even 95 % of the mass. It's 
only the very dilute outer 5 % of the stellar mass that increases 
Rt so much in high mass MS stars. 



2.2 Stellar models 

In our simulations, we use realistic main sequence (MS) 
models to set up the initial stellar str uctures. Models from 
the Geneva stellar evolution group dSchaller et al.l Il992l : 
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iMevnet et"ai] Il994l) have been applied for ZAMS masses 
ranging from 0.8 to 85 Mq , and models bv lCharbonnel et alJ 
l|199fy for masses down to 0.4 M© . For each (initial) stellar 
mass, we had to select one particular model among those 
spanning the main sequence evolutionary track. We chose 
the instant tmodci which divides the MS life in two parts 
with approximately equal collision probabilities. Assuming 
strong gravitational focusing and neglecting any mass loss, 
the collision probability per unit time is dPcoii/dt oc -R*(t), 
so that, 



r mi 



min(tMS.12 Gyr) 



Ri,{t) dt 



(9) 



with t = on the ZAMS. For high-mass stars 20 Mq) 
mass loss by stellar winds is already im portant on the MS 
JSchaller et alll992tlM'evnet et ajigol so that the adopted 
models have real masses lower than their nominal {i.e. 
ZAMS) masses, for instance, the largest star we consider, 
a "85 Mq" model, has an actual mass of only 74.3 Mq. The 
mass-radius relation is shown in Fig. 2] For M ^ 0.4 Mq , 
it is given by the stellar models just discussed. For smaller 
masses, we simply extrapolated a power-law relation from 
the 0.4 Mq and 0.5 Mq points. It appears that this gave 
radii in go od agreement with detailed structure models by 
IChabrier fc Baraffc (1993) who yield ~ 0.12 Rq at 0.1 M© 
and R ~ 0.22 Rq at 0.2 Mq. 

We used models with solar composition {Y — 0.30, Z — 
0.02). A population II metallicity {Y = 0.24, Z = 0.001) 
would introduce significant alterations in the stellar struc- 
tures. Most noticeably, low-Z stars are initially more com- 
pact, with radii smaller by 1 0-40 %, and have larger convec - 
tive cores for > IMq jKippenhahn fc Weiger3ll994^ . 
For high-mass stars, the most important difference i s proba- 
bly t he much weaker mass-loss at lower metallicity (iMaedeil 
Il992fl . We made no attempt to assess the impact of these 
effects on collision outcomes. We hope that they can be par- 
tially scaled out by a proper dimensionless parameterisa- 
tion of the initial conditions and results of the collisions (see 
Sec.|HJ . While the structure of stars less massive than 0.5 Mq 
is very close to that of an n = 1.5 polytrope, more massive 
evolved MS stars do not match any polytropic model. In 
particular, stars with M* ^ 1 Mq are more concentrated 
than 71 = 3 polytropes (see appendix for density profiles). 

The lowest stellar masses considered are 0.1 and 0.2 Mq. 
For such objects we didn't use detailed s tellar structure mod- 
els like those by IChabrier fc Baraffj ^1997^ because they 
rely on a very complex equation of state (EOS) accounting 
for degeneracy and electrostatic effects. Such an EOS was 
not available to us for use in the SPH code at the time we 
embarked on this project. Also, solving this kind of com- 
plicated EOS (for each particle at each time step) is done 
using an iterative scheme and represents a significant com- 
putational burden. Instead, we note that the interior of stars 
with masses lower than ~ 0.4 Mq is nearly completely con- 
vective, so their inte rnal structure is very clos e to that of 
an — 1.5 Dolvtrope illansen fc Kawaleilll994t Chabrier fc 
Baraffe T2000l) . Given the mass and radius, we can build an 
initial polytropic star in hydrostatic equilibrium using the 
EOS for a fully ionised ideal gas. For 0.2 Mq, we have com- 
pared our simple polytropic model with ideal-gas EOS to a 
state-of-the-art stellar structure provided by Isabelle Baraffe 
and found that discrepancies in the density and temperature 



profiles are below 10 % except for the outermost envelope, a 
thin layer which is not represented in the SPH structure. In- 
specting the realistic 0.2 Mq model, we see that only of order 
0.01 % of the stellar mass has temperatures below lO"" K for 
which incomplete molecule dissociation and ionisation may 
be important. Neglecting molecules and partially ionised gas 
may lead to a slight overestimate of the mass loss because 
some of the available kinetic energy has to be used to break 
up molecules and ionise atoms. This is certainly a very small 
effect as the energy required to completely ionise one gram 
of stellar matter of^ola^omgositfoi^^^^^^^^Kippen- 
hahn & Weigert but the kinetic energy at 500 km s 

(a typical contact velocity for a parabolic collision) is of or- 
der 100 times larger. 



2.3 The SPH code 

Smoothed Particle Hydrodynamics is a Lagrangian particle- 
based method that has been widely used to tackle all kinds 
of astrophysical problems, from planetesimal fragmentation 
to cosmological structure formation. For a description of the 
metho d an d of its achievemen ts, we ref er to reviews by|Benz 
"l990l) and lMonaahanl Jl992ft . See also ISteinmetz fc Miillei 



1 19931) for a critical examinati on of the pros and cons of the 
method and iMonaghanI 1^9^ for a presentation of its most 



recent developments. 

We used a versio n of the SP H code that corresponds 
to the description in iBend lll990l) . The kernel function is 
the st andard spline introduced by iMonaghan fc Lattanzid 
1^8^ . This code implements a binar y tree to compute grav- 
itatio nal forces and find neighbours l|Presjll986t iBenz et alJ 
I1990I) . "Bulk" and von Neumann- Richtmyer artificial viscos- 
ity terms are included with a = 2.5 and /3 = 1.0. 

For the stellar matter, we assume the EOS of a com- 
pletely ionised mono-atomic ideal gas with account of the 
radiation pressure: 



2 



T + 



a 

^ 3 
P 



y4 



(10) 

(11) 



where p is the mass density, T the temperature, P the 
total pressure, u the specific internal energy, fi the mean 
molecular weight, a = 7.56 x 10"^*^ J m"^K~* and SH = 
8314JK"^kg"\ The molecular weight of each particle is 
attributed from the initial stellar structure (see next sub- 
section) . It remains constant during the complete SPH sim- 
ulation. In hydrostatic main sequence stars, the radiation 
pressure becomes important fo r masses larger than 5—10 Mq 
ilKippenhahn fc Weiger3ll994^ . 

Release of nuclear energy has been shown to have 
none or very small hy drodynamical infiuence (Rlathis 196i 
iRozvczka et al.lll98Sf) . We thus do not include nuclear re- 
actions in the energy equation. We also neglect radiative 
transport. As long as the gas is optically thick (and the 
bulk of it certainly is during the whole collisional process), 
energy transport by radiation is a diffusion process which 
time scale, for a sun^ik^ta^^^^^^^^^W^eM^Kip- 
penhahn fc Weigert^WQ^^^eMiTHehnho^rTimB^^liis is 
enormously larger that the hydrodynamical time-scales (a 
few tens of hours, at most). For a gas cloud of radius R and 
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mass M, the diffusion time is: 

tdiff ^ with ^ = (j) (12) 

where I is the mean free path of photons. It is connected to 
the opacity k hy I — [np)^^ . Thus, 

,,„^^|.120yrsx(-^)(|^)(^) , 

(13) 

where Kes — 0.04 kg ^m^ is the opacity due to electron scat- 
tering (a lower bound for hi in ionised gas). It follows that 
radiation cooling is negligible even long after the end of the 
collision simulation. 



2.4 Building of initial SPH stellar models 

Building an SPH star from a given stellar structure model is 
not completely straightforward. First, the spatial positions 
of the particles have to be chosen. Then, each particle must 
be given a mass and smoothing length in such a way that 
the total mass is respected and the model's density profile 
p*{R) is well approximated by the SPH interpolate. A sec- 
ond thermodynamical variable (the internal energy u{R), in 
our case), as well as the chemical composition, is also speci- 
fied by the structure model. These quantities determine the 
pressure through the EOS. If the EOS is similar to the one 
used in the stellar structure code, the SPH structure should 
be in hydrostatic equilibrium. 

If all particles are attributed the same mass, their num- 
ber density must closely follow p*(i?), which, unless a huge 
number of particles are used, results in a severe under- 
sampling of the outer regions where the "action" takes place 
during most collisions. On the other hand, using a con- 
stant particle density throughout the star, by placing them 
on a periodic grid for instance, will lead to a very inaccu- 
rate core representation as a small set of very massive par- 
ticles. This could expectedly yield unstable initial models 
and noisy collisional results in case these few heavy particles 
strongly participate to the hydrodynamics of the encounter. 
We thus had to find some compromise between these two 
extreme options. If we neglect particles overlap, the relation 
p^(7i) ~ mpart(-R)npart(-R) holds, with m-part being the iocal 
mass of each SPH particle at distance R from the star's cen- 
tre and Mpart their number density at that position. Thus, 
we decided to impose 

ripart (X and mpart oc (14) 

with a = 0.5 (the two above mentioned extremes correspond 
to Q = 1 and a = respectively). To obtain a i?-dependent 
w-part , we place particles on concentric spheres with variable 
spacing. On each sphere particles are arranged on constant 
"latitude" circles. When this is done, the smoothing lengths 
hi are adjusted until each particle overlaps approximately 
with the same number (~ 40) of neighbours' centres. Fi- 
nally, particles' masses are iteratively adjusted in order to 
bring the SPH interpolate for the density (at the centre of 
particles) closer to the model's p*. This is done by repeating 
the assignments 

m. ^m, (0.3 + 0.7 ^*^^'^ ^ i=l...iVpart 



Mi=1.2Mq M2 = 5.0Mg V"^j=6.75xl0=kms-i 




0.2 0.4 0.6 



d^in/(Ri + R3) 

Figure 5. Fractional losses of mass (top panel) and energy (bot- 
tom) as functions of the distance of closest approach (in Kcplcrian 
approximation) for fixed (Mi, M2, V^^)- Three sets of simulations 
are reported. In the first one (solid lines) we use a relatively low 
number of SPH particles and our standard ( rather large) valu e 
for the binary tree accuracy parameter acc jBenz et al.lll990ft . 
In the second series, we used about four times more particles. In 
the third, we used a lower acc value which gives a more precise 
computation of gravitational forces. The results are nearly coin- 
cident. Tcont is the kinetic energy at contact, i.e. for a separation 
between centres of Ri + R2 ■ 



20 times. As this procedure doesn't conserve the total mass 
M* exactly, all are then slightly rescaled by a uniform 
factor to obtain the required M, . This method is fast and 
effective to give a good match to for the bulk of the stellar 
interior, as testified by the profiles shown in the appendix. 
But, despite the use of lighter particles to represent the gas 
in the stellar envelope, the outermost layers of the star are 
poorly modelled. In particular the SPH realisation fails at 
precisely reproducing the stellar radius. This had to be ex- 
pected in models with a limited number of particles. 

Nonetheless, in grazing collisions, our use of low-mass 
SPH particles to represent the outer parts of a star appar- 
ently leads to a reliable determination of fractional mass 
losses as small as 10~'*-10~^. This claim is grounded on 
diagrams like Fig. |^ which shows the fractional mass loss 
for two sets of simulations, the first one with the "normal" 
(low) resolution and the second one with a number of par- 
ticles about four times larger. The differences are obviously 
very small for all cases but the most distant interactions. 
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An extended coverage of the four dimensional 
(Ml, M2, Koi, b) parameter space requires a huge number of 
collisions to be computed. On the other hand, we don't need 
very accurate results; a relative precision of about 10 % on 
mass and energy loss should be sufficient for our purposes. 
More precise results would not make much sense anyway as 
any application will probably require some sort of inter- or 
extrapolation from our simulation data, to apply it to stars 
with different masses or metallicities, for instance. Thus we 
decided to tune numerical parameters to values that allow 
relatively fast computations while ensuring reasonable ac- 
curacy. This means that we generally used 1000-8000 SPH 
particles for each stellar model (a few collisions have been 
computed with the most massive star having ^^16 000 or 
~32 000 particles). In most simulations the total number of 
particles ranges between 2000 and 10 000. Thus, a collision is 
computed in a few hours to a few days on a run-of-the-mill 
workstation. We use a higher number of particles in high 
mass stars in an attempt to resolve both the high density 
centre that contains most of the mass and the low density 
envelope which is more likely to interact with the other star. 
We also adapt the number of particles of both stars in or- 
der to get spatial and mass resolution not too dissimilar for 
stars of unequal sizes. As an example, for equal mass stars we 
generally use 2000-1-2000 particles for low masses (sj IM0) 
and 4000+4000 particles for high masses. These numbers 
are certainly not impressive by present-day standards but 
already corresponded to considerable computational burden 
given the large number of collision to simulate and the typi- 
cal speed of computers at the time this study was initiated, 
more than five years ago. We now discuss the test com- 
putations we made to ensure these particle numbers were 
sufficient for our aims. 

2.5 Determination of the required resolution 

To determine the minimal desirable number of particles to 
be used in our simulations, we computed the same physi- 
cal collision with various A^'part. Fig. 10] shows the evolution 
of the energies during two typical collisions simulated with 
increasing numbers of particles. In all these cases but one, 
the number of particles in the small star is ~1000 while the 
large star consist of ^2000 to ~32 000 particles. The first 
collision is a high- velocity quasi-hyperbolic "fly-by". In the 
second example, the stars are left bound to each other after 
the first periastron passage. A second collision ensues that 
leads to a rapid merging of the small star into the centre of 
the big one. In the "fly-by" case, the energy evolution curves 
for the various A^'part values are very similar to each other 
as soon as more than 2000 particles are used to represent 
the large star. Contrariwise, in the "merger" case, the time 
between the two successive periastron passages exhibits a 
strong dependency on A^part. Not only does it not converge 
to some "real" value as the resolution is increased but the 
opposite occurs! This intriguing behaviour casts important 
doubts on the ability of this version of the SPH code to follow 
reliably the formation and evolution of tidal binaries. How- 
ever we note that a simulation with 32000 -|- 2000 particles 
(the last one in the right panel of Fig|SJ gives nearly exactly 
the same energy evolution as another one with 16000 -I- 1000 
particles (the fourth one). This is a hint at the importance 
of a good resolution in both stars and not only the larger. 



Convergence in the results is attained only if we increase 
both resolutions. It is not clear to us why this is so but it is 
obviously connected to the poorly resolved envelope struc- 
ture which determines how much orbital energy is dissipated 
at first passage. This turns out to have very little implica- 
tion for the amount of mass loss, the only quantity we want 
to determine in case of a merger. As shown on Fig. |K| it 
is only for grazing fly-bys that one get significant relative 
discrepancy between different resolutions, as long as energy 
and mass loss are concerned. This also is due to inaccura- 
cies in the SPH realisation of the outer layers. But given 
the very small absolute values of these losses, this can only 
have very small impact when SPH results are used to imple- 
ment collisional effects in simulations of high- velocity stellar 
systems. 

In any case, in our work, we are mostly interested in the 
final outcome of collisions in terms of global quantities, like 
the mass and energy loss. The dependency of these quanti- 
ties on A'part turns out to be very weak. The fractional mass 
and energy losses and the fractional deviation from the Ke- 
plerian defiection angle typically vary by less than 20 % over 
the whole range of particle numbers used in these tests (see 
appendix for diagrams). While the lowest A^part produce re- 
sults that are somewhat off, as compared to higher resolution 
runs, 1000-1-2000 particles seem to be already sufficient. 

2.6 Starting, ending and analysis of a collision 
simulation 

Any collision to be computed requires the specification of 
both stellar models and of the relative velocity at infinity 
Vj.^ and the impact parameter b. We neglect any finite- 
size effect until the separation between the stars' centres 
is 3(-Ri -I- -R2). Hence we analytically advance the stars to 
this situation on hyperbolic trajectories. At this point, we 
start the SPH simulation and set t = 0. 

The computation stops at tend ^ 20rdyn with Tdyn = 
^y{Ri + i?2)V(G(Mi -I- M2)). If, at this time, two surviving 
stars are present with separation less than 3(-Ri -|--R2)ini or if 
the amount of gas with uncertain fate (see next paragraph) 
exceeds 1 % of the total mass, the simulation is integrated 
further (by setting tend <— 2tond) until it passes these tests 
or some maximal integration time is reached. In practise, no 
collision required integration past t = 25OO\/R0"^/(GM0). 
Unfortunately, these simple termination prescriptions are 
probably inadequate when a bound binary forms after first 
periastron passage. They can indeed force integration for 
many orbital periods although we expect the SPH scheme 
(at least when used with our set of numerical parameters) to 
lose reliability in that regime whose outcome is very likely to 
be a merger (see Sec. 13. H . A wiser approach would have been 
to identify such encounters, stop computation after first pas- 
sage at periastron and, if needed, to rely on other theoretical 
considerations to assess the outcome. 

One monitors the energy (non-)conservation with Se = 

-Eond — -Einil / EnoTui with -Enorm ~ -Ekhi + -E-bind, whcrC -Eini, 

Ecnd are the initial and final total energies, the initial 
kinetic energy at infinity (orbital energy) and -Ebind the sum 
of the binding energies of the stars (positive definite). Using 
the total energy, i?ini = -Ekin ~ -Ebind for normalisation isn't 
appropriate because it may be very close to zero, leading to 
misleadingly large values oi Se- Enoim gives a natural energy 
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Figure 6. Evolution of the system energies during two collisions. Each collision has been computed with different resolutions (3000 to 
34000 particles). Code units are used here. The quantities describing the collisions are specified on top of the diagrams. See text for 
further comments. The left column is a "fly-by" encounter. The right column corresponds to a merger (see Fig. 171. 



scale for the problem. Using -Bbind for normalisation doesn't 
change much the 5e statistics. The worst non-conservation 
amongst all simulations is Se — 0.06; all but 15 simulations 
have Se < 0.01; 66% of all runs have Se < 0.001. 

The SPH code yields quantities describing every particle 
at the end of the computation, i.e. their positions, velocities, 
internal energies. . . This raw "microscopic" data has to be 
analysed to provide a useful description of the outcome of 
the collision in terms of "macroscopic" quantities, i.e. prop- 
erties of outgoing star(s) (if any). Namely, we want to know 
how many stars survive (0, 1 or 2) and what their masses, 
positions and velocities at the end of the computation are. 
This, and any other aspect of the structure of the star(s) (see 
Sec. 13.311 and of the ejected gas, can be easily determined if 
we manage to build a list stating to which star each particle 
belongs or whether it is unbound to any surviving star. This 
data is provided by an analysis algorithm which proceeds in 
two steps: 

(i) A first guess attribution is realised by a code which 
tries to identify, through density and proximity criteria, zero, 
one or two concentrated lumps of particl es. To this end, we 
use th e freely available HOP algorithm bv lEisenstein fc Hud 
I^Q^)- These groups are regarded as first approximations 
of stars to be refined in the second stage of the method. 



(ii) We then iteratively cycle through all particles to com- 
pute the energies of each one relative to both stars ( "A" and 
"B"). For iteration k, the energy of particle i relative to star 
A as identified at the previous iteration {"Ak-i") reads 

aI 1 / ^ \ ^ 

-Ei 1^ = Ui + -mi l^Vi — VAk_ij 

-Gm. E (15) 

In this formula, Ui, m;, Vi and Xi are the internal (thermal) 
energy, mass, velocity and position of particle i. VAt,_i is the 
velocity of "star" A^-i, i.e. 

If either E^\^. or Ef\^ is negative, the particle is ascribed 
to the "star" relative to which its energy is the most neg- 
ative. Particles with positive energies relative to both stars 
but with negative total energy in the collision centre-of-mass 
reference frame (CMRF) are tagged as "doubtful". At the 
end of the iteration we thus get new sets of particles making 
up "stars" Afe and B^. We go on iterating until no modifi- 
cation in the composition of these sets occurs anymore. 
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This procedure deserves a few comments. 

• The energy criterion may fail to predict the correct at- 
tributions. For instance, a particle with high velocity toward 
a given star may happen to have positive energy relative 
to this star even if it will impact it and thus merge into 
it. Furthermore, even without resorting to hydrodynamical 
processes, we learn from studies of the gravitational 3-body 
problem that the eventual fate of a particle submitted to the 
gravitational forces of two massive bodies can not generally 
be predicted just through energy consideration. However, if 
we carry on the SPH integration to a physical time large 
enough for the stars to have moved away from each other to 
a large separation and/or for the large amplitude hydrody- 
namical processes to have ceased, we expect the final SPH 
configuration to be essentially free of such problematical par- 
ticles and the energy criterion to be reliable. 

• "Doubtful" particles generally lie between the two stars 
so that they gain negative gravitational contributions to 
their total energy from both potential wells even though 
they are not bound to any one star. In such cases, their 
number should decrease as the distance between the stars 
increases. Another situation that can leave a relatively im- 
portant doubtful mass fraction (i.e. > 1 % of the total mass) 
occur in high velocity head-on collisions that result in an 
expending gas cloud. Its central part, lying in the potential 
well of the surrounding gas has negative total energy but 
nevertheless expands to infinity. Although these cases seem 
to have genuine physically interpretations, there are situa- 
tions where a high Mdoubt is indicative of some error in the 
analysis. One such case consist of a close tidal binary being 
erroneously identified as a single particle group in the first 
step. The iterative steps then progressively interpret one of 
the stars as a group of doubtful particles while retaining 
only the other lump as a "real" star. 

• This last example illustrates how critical the first attri- 
bution stage is. Its failure to detect independent stars can- 
not be recovered by the iterative process! There is probably 
room for improvement in this part of our analysis procedure 
and the use of HDP, an algorithm aimed at finding structures 
in large cosmological simulations, is arguably an inefficient 
overkill. A simple-minded approach that first divides parti- 
cles in the same two groups that built up the pre-encounter 
stars proves to allow convergence to real stars in cases that 
confuse HOP. To account for mergers, when the distance be- 
tween the centres of the two groups is much smaller than 
some typical size, we can ascribe all particles to a single 
group to be then "eroded" down to the bound star (if it 
exists) by the iterative energy test. 

All simulations were first analysed using HOP to produce 
the initial particle attributions. The results of the iterative 
procedure were then visually inspected by plotting log p ver- 
sus spatial coordinate x for all SPH particles and using dif- 
ferent colours to code the attributions. Errors are immedi- 
ately spotted in such a diagram allowing one to integrate 
the simulation for a longer time if the separation between 
"stars" (density peaks) is deemed too small or switching to 
the just-mentioned simple-minded scheme for initial attri- 
butions in the few cases HOP clearly made a wrong guess. 

In the vast majority of simulations, we only run the 
analysis software just described on the final SPH file. As 
mentioned above, if, for that configuration. A/doubt exceeds 



some fraction of the total mass (1 %) or wrong attribution is 
seen, we compute the interaction for a longer physical time. 
When the integration is deemed over and the properties and 
kinematics of the surviving star(s) have been determined, we 
assume that the stars' masses have reached constant values 
and that the subsequent orbital evolution is purely Keplerian 
again. This allows to compute Scoii as an asymptotic value. 
The physical time tond over which the SPH simulation is 
computed has thus to be long enough for the strong hydro- 
dynamical regime to be over. On the other hand, choosing 
too large a value for tend, is not only computationally ex- 
pensive but could result in inaccurate results due to the ac- 
cumulation of small numerical errors. Hence, it is of interest 
to analyse a few typical collisions at a number of increas- 
ing times during the SPH computation to test whether the 
outcome quantities have reached steady values and whether 
these values show sign of numerical drift at large t. Fig. |7| 
is an example of such computations. The plot of the trajec- 
tories (panel (a)) testifies that, in most cases, the analysis 
procedure identifies the stars correctly, even during close in- 
teraction. The curves for the evolution of predicted mass 
and energy losses show abrupt increases at periastron pas- 
sages and stay nearly constant quickly after the last close en- 
counter (leading to a merger) is over. Although the analyse 
software gets confused when the stars penetrate each other, 
this is of no practical concern because it is only a transi- 
tory situation. For fly-bys (including the case the small star 
emerges as an unbound, expanding cloud), we integrate until 
the stars are again very well separated; when stars capture 
each other, the analysis is only done after a merged object 
has formed or when the stars, forming a binary, do not over- 
lap. We conclude that the way we terminate SPH collisions 
and analyse their results is sound. 

2.7 Building a comprehensive table of collisions 

This study was first embarked on as a sub-project. It is part 
of a work aimed at simulating the stellar dynamical evolu- 
tion of dense galactic nuclei hosting black holes. To this end, 
a new Monte Carlo code for cluster dynamics has been devel- 
oped [Frcitag & Bcnz 2001b. 200^) . In order to incorporate 
the effects of stellar collisions with a high level of realism 
into it, we decided to compute a large number of SPH sim- 
ulations spanning all the relevant values of the initial con- 
ditions. Our hope was then to extract fitting formulae from 
this database of results to get an efficient description of the 
outcome of any arbitrary collision that could occur during 
a cluster simulation run.* Figuring out such mathematical 
descriptions proved too difficult and we have resorted to an 
interpolation procedure. This will be explained in Sec. 13.31 
Contrary to globular clusters where all collisions are es- 
sentially parabolic due to the velocity dispersion being much 
lower than the escape velocity from a stellar surface, galac- 
tic nuclei may have deep potential wells, or even harbour 

* A collision requires a few hours to a few days of CPU time to be 
simulated by the SPH code on a standard workstation and some 
simulated high density nuclei experience many thousands of these 
events during a run spanning a physical duration of lO'^'^ years. It 
is consequently clearly impossible to switch to on-the-fly SPH in- 
tegrations when collisions are detected in the cluster simulations! 
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Figure 7. Collision between stars with Ah = IMq, M2 = SM©, V;^ = 0.07V; = 43.71cms-\ dmi,i/(-Ri + R2) = 0.39. Panel (a): 
Trajectories of colliding stars, as identified by the algorithm used to analyse the outcome of collisions. This encounter leads to the 
formation of a binary which coalesce after two orbits. The analysis algorithm is unable to tell one star from the other during the final 
merging. Panel (b), top: Evolution of the separation between both stars. Panel (b), bottom: Evolution of the amount of gas unbound 
to any star, with positive or negative energy in the centre-of-mass reference frame (5Mioag and <5J\/doubt respectively). This diagram 
illustrates how the mass loss increases abruptly at each periastron passage and reaches a steady value after complete merging. The 
amount of gas with doubtful fate also gets quickly to a vanishingly small value. This ensures that the interaction has been integrated for 
a sufficiently long time. 



massive black holes and thus force some of their inhabiting 
stars to collide on high-velocity hyperbolic trajectories. For 
instance, at the centre of the Milky Way, "SO" stars are 
on orbits with pericentre velocities of up to ~ 12 000 km s~^ 
ijGhez et alJl200!Tl: ISchodel et ai]l20n.1) and even higher val- 
ues will probably show up in future higher resolution obser- 
vations reaching closer to the ~ 3 — 4 x 10^ Mq black hole 
Sgr A* . Hence, we cannot restrict ourselves to collisions with 
V^'^ ~ but have to go up to a few thousands of kms~^. 

Moreover, the population in galactic nuclei does not 
consist of old stars all born at the same time but may include 
MS stars with an extended range of masses. High mass stars 
are particularly important in the first phases of the system 
evolution: relaxation-induced mass segregation may quickly 
concentrate them in the high density central regions where, 
having large cross sections, they get relatively high collision 
probabilities, despite their overall scarcity and their short 
MS lifetimes. Consequently, we have also to span a large 
range of initial masses, extending far beyond the ~ 1 M0 
tun-off mass that would be sufficient for a study of colli- 
sions in present-day Galactic globular clusters. 

Finally, to further extend the domain in parameter 
space to be explored, we note that stars of different masses 
have very different internal structures (see density profiles in 
appendix) so we cannot hope to scale out the absolute mass 
from the collision process. For instance, we would expect 
the (dimensionless) results to depend only on the mass ratio 
M1/M2 only if stellar structures were homologous to each 
other and a power-law M-R relation was obeyed. As this is 



not true, we have to consider the masses of the in-coming 
stars as two independent variables. 

Summing it up, we have to deal with a fully 4- 
dimensional parameter space in which we have sampled a 
domain which is more or less the following: 

• Stellar masses from 0.1 to 74.3 M0 (the latter value 
corresponds to a ZAMS mass of 85 M0). 

• Relative velocities in the range /V, ~ 0.03-30. 

• Impact parameters corresponding to dmin/(_Ri +-R2) = 
0-0.9. 

A mere 10 points resolution for each parameter already 
turns to a total of 10 000 collisions to be computed, a num- 
ber clearly beyond what can be managed "by hand". This 
high number grounded our decision to neglect other, "second 
order" parameters such as metallicity, rotation or evolution- 
ary stage along the MS track. A complete software system, 
consisting of many UNIX shell scripts has been developed 
to run these SPH simulation in a (nearly) completely auto- 
matic way. The system looks through a table for collision 
simulations that have not yet been computed to their end 
and makes them run on idle workstations available through 
the local computer network. The system interrupts a simula- 
tion job when the computer on which it's running ceases to 
be "available" (basically during daytime) and calls the anal- 
ysis software when a run is over. If no further integration is 
required, the results are added to an output table. Supervis- 
ing this automatic system is not as painless as it may sound: 
due to the number of simulations that run concurrently (10 
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14 M. Freitag and W. Benz 



to 50, typically) , "exceptional" problems mainly originating 
from malfunctions in the local network occur nearly every 
day and have to be fixed manually. All in all, obtaining a 
system reasonably crash-proof revealed itself to be unex- 
pectedly difficult. This paper reports on the results of the 
~ 14 000 simulations we managed to compute with this ap- 
proach. On Fig.|H| we attempt to show the initial conditions 
for all simulations. 



2.8 Formation of binaries through tidal 
interactions 

Even when the periastron distance is larger than the sum of 
the stellar radii, close encounters at low relative velocity can 
rise tides in the interacting stars and lead to the formation 
of a b ound binary. As already pointed out bv lFabian et alJ 
1^7^, in globular clusters, the cross section for such tidal 
captures is a factor 1-2 times as large as for collisions (as- 
suming a typical relative velocity of lOkms^^). Determin- 
ing through SPH simulations the critical impact parameter 
for tidal captures in (quasi-)parabolic, non-touching encoun- 
ters is a demanding task, requiring high resolution of the 
stellar envelopes where tides transfer energy from the or- 
bital motion to stellar oscillations. This phenomenon is not 
treated in this paper because, in typical galactic nuclei, the 
relative velocities are in excess of 50kms~^, a regime where 
tidal binaries can form only for very close encounters, re- 
quiring contact interaction in most cases, with the possible 
excegtion^^fis^on^en^^te^Jo^^nas^^tf^Lee & Nel- 
son 119881: iKim fc Leelll999^ .' Hence, we restricted ourselves 
to the range dmin < {Ri + Ri)- 



3 RESULTS 

3.1 Overall survey of the results 

Trying to get a complete coverage of collision parameter 
space implies a huge volume of simulation results. The dif- 
ficulty of our approach is to extract useful information in 
manageable form out of these data. As the database was 
nearing completion, we looked for mathematical relations 
between various input and output quantities. Due to the de- 
terministic nature of collisions, many strong correlations are 
clearly visible but finding fitting formulae for them eluded 
us. The basic difficulty stems from dimensionality of the ini- 
tial parameter space which seems to be genuinely 4D. 

Here we do not show the results from specific collision 
simulations nor discuss the physical mechanism at play dur- 
ing them, as this has been done extensively i n previous works 
jBenz fc Hillslll987l Il992l: iLai et aljligga iLombardi et alJ 
[^9^) ~^or the interested reader a few specific simulations 
are presented in the appendix. What concerns us here is a 
description of the simulation database as a whole. 

The simplest, most qualitative, description of the colli- 
sional outcome is the number of outgoing star(s). For given 
initial masses, we can plot a 2D diagram indicating this num- 
ber for all collision simulations performed, as a function of 
the impact parameter and the relative velocity (Fig. [nj. 

Before we comment on that figure, some explanations 
are called for. V^o'liact is an approximate value of the relative 



velocity at "half-mass contact" . It is defined through 

K'olLct = \J iK'sf + {vi'-'Y ■ (16) 

It should be noticed that such a "deep" contact does not 
occur during encounter with large impact parameters; this 
value only serves as a convenient parameterisation that al- 
lows to map the [0, oo[ V,!^ range onto [0, 1[. In these plots, 
each dot represent one SPH simulation. Green dots are colli- 
sions survived by both stars (although significant amount of 
mass loss may have occurred). Blue dots indicate that there 
is only one star left at the end of the encounter. Orange dots 
stand for tidal binaries and red dots for complete disruption 
of both stars. One sees that, for this half-mass parametrisa- 
tion of the initial conditions, the borders between these vari- 
ous regimes are primarily set by the mass ration q — Mi /M2 , 
quite independently of the actual masses. Unfortunately, as 
will be stressed below, this appears generally not to hold for 
more quantitative results. 

These diagrams provide a division of the collisions into 
a few different regimes. Most of the (dmin, Kci) plane is oc- 
cupied by "fiy-bys", i.e. encounters from which two unbound 
stars escape. In some cases, this domain extends to dmin ~ 
like a small wedge between the merger regime (lower veloc- 
ities) and disruptions (higher velocities). It is thus possible 
for a small star to pass right through the centre of a larger 
one and not being disrupted. We detected about 250 such 
cases in our survey, all with g between 0.04 and 0.25 and 
Ml (small star) between 0.7 and 2Mq. Moreover, in about 
one third of these simulations (with 1/ < 1.7), the small star 
gains mass during the interaction while the larger star al- 
ways suffers from important mass loss. It seems even possible 
that in some collisions, the small star, acting like a bullet, 
shatters its target but remains nearly intact. Similar out- 
comes were obtained by BH92 for n = 1.5 polytropes with 
q — 0.2. LRS93 did not find any head-on collision with a 
surviving small star. As pointed out by these author, such 
discrepancies -Eis well as other differences between our re- 
sults and published data, see Sec. 13.21 - probably originate in 
the fact that different stellar models have been used. The 
ratio of stellar central densities is likely to be a key parame- 
ter in allowing such "fly-through" collisions. In all the cases 
identified by us, this ratio exceeds 6. However, the astrophys- 
ical importance of this phenomenon is low because, at large 
relative velocities, collisions with small dmin are unlikely as 
gravitational focusing is quenched. 

Mergers or bound binaries are formed during encoun- 
ters with low relative velocities and impact parameter below 
some critical Amerg- This value depends on the relative ve- 
locity and the masses (mostly through the mass ratio). It 
is apparent as the transition between green and orange or 
blue dots on Fig.|n| It is generally larger than R^^^ + R2^^ for 
f < 0.6 and smaller at larger velocities. An ad-hoc analyt- 
ical parametrisation of A merg as a function o f Mi, M2 and 
u will be published soon jFreitag et alj|2005bh . Remarkably, 
the maximum velocity for a head-on collision that still leads 
to merger is i/ ~ 1.7 — 2.1, quite independently of the stel- 
lar models. The border between this region and the "fiy-by" 
regime at higher dmin is also rather well defined if half-mass 
variables are used. 

All binaries formed in our simulations will presumably 
merge into single stars after a few orbits. The reason for this 
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Figure 9. Number of stars surviving collisions. Red dots are for collisions leading to complete disruption, blue dots for cases with one 
surviving star, orange dots for tidal binaries (very likely to eventually merge) and green dots for encounters writh two surviving stars. 
.See text for further explanations and comments. 
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Figure 10. Fractional mass loss for collisions between stars 
of masses 1 and 5 Mq with a relative velocity at infinity of 
43.7kms~^. We indicate the mass loss for each successive perias- 
tron passage. The line with crosses (green in the colour on-line 
version) shows the mass loss after the first passage, the line with 
open squares (cyan) the mass loss after the second passage and 
so on until the stars have merged. The total mass loss is shown 
by the line with round dots (black). For the most central col- 
lisions (rfmin < 0.425), the evolution was integrated up to final 
merging. This was not possible for more distant encounter due to 
the important increase in computing time this would require and 
the loss of numerical precision to be expected in such long SPH 
integrations. 

is that, at each periastron passage with dmin < (^i + ^2), 
some orbital energy is converted into heat by shocks and the 
stellar radii expand so that at next periastron passage the 
hydrodyna mical interaction is stronger and more energy is 
dissipated llBenz et alJIlQSsI) . Hence, the fate of these bina- 
ries is not as complex an issue as the long-term orbital evo- 
lution of systems f ormed through more distant encounters 
jMardlinelll995allJ) . Thus, the border between the regions 
of merging and binary formation probably results from the 
criteria we use to stop the SPH computations and has no 
strong physical meaning. If it were possible to integrate the 
evolution for many orbital periods, there is little doubt that 
any binary would eventually merge. Fig. 1101 illustrates this 
point. To produce this diagram we computed a set of col- 
lisions with for given stellar models and a 
fixed value of which is sufficiently low that every colli- 
sion leads either to direct merger or binary formation. Unlike 
the bulk of our simulations for which we analysed only the 
"final" state, here we report the mass loss after each succes- 
sive periastron passage. Obviously, as rfmin is increased, the 
number of orbits preceding the final coalescence get higher 
and higher, as does the orbital period. Consequently, the 
required CPU time grows up to unacceptable values. A no- 
ticeable feature of Fig. lini is that all the collisions apparently 
converge to nearly the same total mass loss at merging. The 
reason for this behaviour is unknown to us. 



Apart from the low velocity merging zone, another re- 
gion with one surviving star is present in the diagrams of 
Fig. |U] This second zone is more or less confined between 
cases where stars are completely destroyed (for lower impact 
parameters) or both survive (for higher impact parameter). 
This "one-star band", which does not show up when the 
two stars are (nearly) identical, is populated by collisions 
during which the small impactor is disrupted without being 
accreted into the large star. In such high velocity collisions, 
the small star accumulates so much thermal energy as it flies 
through the massive one, that it turns into an unbound, ex- 
panding gas cloud. 

The most spectacular collisions are those that lead to 
complete disruption of both stars. However, to achieve this 
result, we note that both a high relative velocity and a small 
impact parameter are required, a combination made unlikely 
by the absence of gravitational focusing at such high veloc- 
ities so it is clear that neither mergers nor complete disrup- 
tions are likely outcome in galactic nuclei, as confirmed by 
Monte Carlo simulations iPreitag fc Benzl2002tlFreitag et alJ 

liooii). 

3.2 Comparison with literature 

In this section, we perform critical comparisons between 
our results and data and methods previously published (see 
Sec.O- 

The first attempt at quantitatively predictin g the o ut- 
come of off-axis stellar collisions was presented bv lSSGsj As 
it is both elegant and simple (but also very approximate), 
we implemented it for comparison purposes. This allowed us 
to apply it to the same stellar models that we used in SPH 
computations. With no particular optimisation or numer- 
ical tricks, this algorithm computes the results of 50 stel- 
lar collisions in less than 3 seconds on a standard work- 
station! In comparison, a typical SPH run takes about one 
day of CPU time. In our version of this method , which is 
nearly identical to that of iMurphv et al] il99ll) . we con- 
sider that the stars encounter on straight line trajectories 
with an impact parameter (distance between parallel tra- 
jectories) set to dmiii (Eq. 2J and a relative velocity equal 
to Koi = Knax = KS'(b/rfmin) (scc Scc. O for the defini- 
tions of these quantities). We then proceed by dividing both 
colliding stars into small sticks of square cross-section that 
are parallel to the trajectory. The result of the collision, in 
terms of mass and energy loss, is computed by considering 
completely inelastic (i.e. "sticking") collisions between one 
mass stick from each star in the overlapping cross-section. 
Stick i of star 1 collides with stick j of star 2 if they have 
coincident position in the plane perpendicular to the rec- 
tilinear trajectories. No energy or momentum exchange is 
taken into account between stick i and other mass elements 
from its "parent" or the other star. We further assume that 
all kinetic energy to be dissipated to merge i and j is con- 
verted into thermal energy to be shared between these two 
elements and that there is no heat exchange between them 
so that the thermal energy is given to i and j in proportion 
to their pre-collision kinetic energies in the collision centre- 
of-mass frame. Finally, the condition for mass element i to 
be liberated is that its acquired specific energy is larger than 
the initial escape velocity of its parent star, . As demon- 
strated bv lMurphv et al.1 (Il99ll) . this results in the following 
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Figure 11. CoUisional fractional mass loss. We compare some of our simulations (dots and solid line spline interpolation) with results 
from the literature (see text). To obtain the "Spitzer &; Saslaw 66" {SS66) curves we applied the method of these authors to our stellar 
models. Panels (a) and (b): for such small mass stars, the structure is quasi-identical to a n = 1.5 polytrope. This is why we get a good 
agreement with BH87 and BH92 but big discrepancies with formulae from LRS93 and R99, as these authors use more concentrated 
n = 3 structures. Note that the SS66 prescription gives very satisfying prediction for off-axis encounters! Panel (b) This is a case with 
relatively good agreement with published formulae. Still, the predictions from R99 and LSR93 are 2-3 times larger than our mass losses. 
The agreement with SS66 is excellent as soon as rfmin/(^i -|- ^2) > 0.15. Panel (d): Here, the best agreement is obtain with the R99 
formula, despite the velocity being a bit lower than the range explored for this work. SS66 gives satisfactory results, but not LRS93. The 
reason for this discrepancy is unknown. See caption of Fig. 1141 for the probable explanation of the bump at dmin/ (Ri + R2) — 0.2 in our 
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simple escape condition for element i of star 1: 

> ^ (17) 
Ami + Aruj VtcI 

where Arriij are the masses of sticks i and j. For a given col- 



lision, we iterate this procedure a few times with increasing 
resolution (decreasing the cross-section of the sticks) until 
the result converges to some prescribed precision level. As 
can be judged from this description, the number and im- 
portance of simplifications in this approach is impressive. 
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It is thus difficult to figure out the regime(s) in which we 
expect them to hold true. The assumptions on rectilinear 
motion, the use of Vlnax, and the escape criterion leave lit- 
tle hope that sensible results can be obtained either for 
low velocity encounters, or for nearly head-on collisions or 
for cases where high fractional mass-loss is expected (high 
Kei but small dmin)- In an attempt to get better prediction 
at low impact parameters, we im plemented the following 
trick, inspired by ISanderj l)l970a|) . For each star, a "core 
radius" is defined; it it the radius enclosing 1/4 of the to- 
tal mass. An "effective" transverse distance is used instead 

of dmin, dcff = min(dniin, -Rcoro.l + Rcorc,2)- ^cS IS UScd tO 

determine the overlapping sections of the stars and to set 
the effective relative velocity during the collision, through 
Vli = {V.fif + 2G(A/i -f Af2)/doff. This recipe is admit- 
tedly quite arbitrary and, if SS66-like treatment of collision 
is to be used in stellar dynamical simulations, one should 
experiment with other similar prescriptions to find the most 
satisfying one. 

All the other literature results included in our compar- 
ison were obtained through SPH simulations. The pioneers 
in this field were Benz fc Hilla a987..1992i') . They did not try 
to describe their results with fitting formulae, so we can only 
compare t heir simulations to cases with very similar initial 
conditions. iLai et all il993h performed a more extended nu- 
merical survey from which they devised a general empirical 
mathematical description to represent the fractional mass 
loss as well as the critical dmin for merger/binary formation. 
Although it is already clear from the figures of their paper 
that this all-encompassing fit does not provide a very pre- 
cise adjustment of their mass-loss results, we use it anyway 
for our comparison. This permits an assessment of the util- 
ity of such formulae as an interpolation tool. To the best of 
our knowledge these formulae have never been adopted to 
incorporate the effect of collisions in stellar dynamics sim- 
ulations. For his study of the coUisional evol ution of a star 
cluster bound to a supermassive black hole, iRauchI (^9^ 
derived another set of fitting formulae from a set of collision 
simulations performed by Melvyn Davies. Individual results 
from these simulations are not published but it's worth men- 
tioning that Rauch-Davies' formulae give not only the mass 
loss but also the energy loss and the (non-Keplerian) angle 
of deviation for the trajectories. 

These comparisons are motivated by two complemen- 
tary goals: 

(i) To test our results. Although the SPH code has be 
throughly tested in the past, we had to develop new tools 
for the present work. For instance, we developed the program 
to compute initial conditions and stellar structures and the 
one that carries out the analysis of the stellar outcome at 
the end of the simulation. To perform this check, we have 
to choose, in our runs and in the literature, cases that have 
initial conditions and stellar structures agreeing as closely 
as possible with each other. 

(ii) To assess whether already published results, which 
covered only a limited region in the parameter space, still 
yield meaningful results when extended beyond this zone. 
We thus dare to compare some of our results with data ob- 
tained using quite different stellar models or with prediction 
of formulae that we apply outside the parameter domain for 
which they were established. Such confrontations should cer- 



tainly not be seen as a way to cast doubt on those published 
results but as an a posteriori motivation for our own work. 

All our comparisons focus on the fractional mass loss. 
This quantity is presumably the most important for inclu- 
sion of the effect of star-star collisions in stellar dynamics 
models and it is given in all previous works. In Sec. 13.31 we 
explain that, in a general case, the description of the out- 
come of a collision requires at least 4 quantities. 

In Fig. 111! we show some selected cases for which we 
expect a good agreement with the literature results. There 
are however some exceptions that we explain in the cap- 
tion of this figure. In Fig |12l more extreme comparisons are 
made. With Fig. [HI we concentrate on comparisons with 
simulation results of LRS93. 

After inspection of these plots and many others not 
shown here, the following comments can be made: 

(i) When comparing some of our simulations to other in- 
dividual computations with very similar initial parameters, 
a comforting, if not surprising, agreement emerges. This is 
particularly true of results from BH87 and BH92 ^. The sit- 
uation with LRS93 (Fig. [HJ is more complicated and we 
discuss it in detail below. 

(ii) The initial stellar structure plays a central role in de- 
termining the results. But this strong dependency may prob- 
ably be compensated to a large amount by some "clever" 
parameterisation of the initial conditions (see below). 

(iii) Fitting formulae can not be used as extrapolation 
tools. This means not only that we cannot trust them when 
applied to larger or smaller velocity, masses or impact pa- 
rameter values than the ones the have been forged for, but 
also that they will fail at predicting outcomes for other stel- 
lar models. 

(iv) Predictions from LRS93 and R99 formulae are gen- 
erally quite different, even when applied to the parameter 
domain in which they should both be relevant. This may 
be due to variations in the stellar structure (the M-R re- 
lation) and/or amplified from small differences at the SPH 
level by the fitting procedures themselves. This is another 
indications that such formulae should be use with extreme 
caution. 

(v) An unexpected result from these confrontations is 
that the best match at dmin/(J?i +^2) > 0.15 and V^'^ ^ 1 is 
obtained with the SS66 method, which incorporates nearly 
no real physics! Furthermore, some of the crudest assump- 
tions it relies on, which are certainly to be blamed for its 
breakdown at low impact parameter may probably be im- 
proved on. An exploration of the real potentialities of this 
simple approach would be an interesting follow-up of the 
present study, mainly because it reduces stellar collisions 
to very simple considerations about momentum and en- 
ergy conservation and could thus be a useful guide toward 
a deeper insight into these processes. Once again, this un- 
expectedly good agreement strongly hints toward the cen- 
tral importance of the stellar structure in collision simula- 

^ BH87 made use of an earlier, much simpler version of our SPH 
code. The smoothing length had a unique, non-evolving value, an 
exponential kernel was used and the gravity was computed by 
direct summation. The code used by BH92 included essentially 
the same features as ours but all particles had the same mass (as 
in BH87). 
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tions. We should add that the SS66 approach also appar- 
ently breaks down for very high velocities V^'^ > 10 where it 
yields too small a mass loss as compared to our simulations. 
It is interesting to note that the parameter domain for which 
SS66 gives very good results is well suited for collisions oc- 
curring in dense galactic nuclei. It may thus be that this 
recipe, when complemented with some prescription describ- 
ing the domain of complete disruption, can be made into a 
useful ingredient for the study of such systems. 

Let's now focus on to comparison with results of LRS93, 
illustrated by Fig. 1131 This work is of special importance 
as it constitutes the only survey of some breadth, also in- 
cluding high-velocity encounters, published so far. These au- 
thors used Eddington models with n = 3 polytropic den- 
sity profiles and assumed _R, oc M,'^. Eddington models 
have a constant (3 = Pgas/^tot; they can be parametrised 
by a = 7.89(1 - P^^^p-^. LRS93 further assume q = 
M^/My = a</a> ^ 1 where subscripts < and > indicate 
the more and less massive star in the encounter, respec- 
tively. LRS93 have parametrised their results through a set 
of formulae that give the fractional mass loss as a function 

of q, a>, Woo := y,S(GM>/i?>)-l/2 dn,in/{R< + R>). 

When comparing our results to this parametrisation, we set 
/?> := 2Ey /Wy where Ey is the total energy of the massive 
star (thermal plus gravitational) and Wy the gravitational 
contribution. This relation is exact for Eddington models 
and is used here to define some "effective" /3 parameter. /3 
is very close to 1 for M* < 10 Mq , leading to small a values. 
Hence, most LRS93 results (with a> = 10, panels a and c 
of Fig. are adapted to My > 10 Mq . This probably ex- 
plains why LRS93 get considerably more mass loss than us; 
their stellar model have little binding energy compared to 
our realistic MS stars. Indeed, the best agreement is reached 
with the few a> = 1 models, see Fig. 113b . d. Also, we stress 
again that n = 3 polytropes do not represent in a satis- 
factory way the mass distribution of any (evolved) MS star 
except, maybe, around M, = 0.9 (see appendix). 

We now turn to an examination of the impact of the 
stellar models on collision results. In Fig. 1141 we compare 
two sets of simulations. In both series, we computed colli- 
sions between stars of masses 0.5 and 2.0 Mq for two differ- 
ent relative velocities and a sequence of impact parameters. 
In the first set, we used realistic stellar models, while in the 
second series, the small star is represented as a n = 1.5 poly- 
trope (which is a very good approximation) and the large 
one as a n = 3 polytrope (a poorer model) . Panel (a) in this 
diagram strongly confirms point [(ii)| of the previous enumer- 
ation. Except for head-on encounters, the polytropic models 
systematically overestimate the mass-loss by factors as large 
as 5! This seems to strongly justify our use of realistic stellar 
structure instead of the traditional polytropes but panel (b) 
slightly weakens this statement. There we use the half-mass 
radii to normalise c?min. This simple change of parameter, 
scales out the discrepancy to a large amount. Only for large 
dmin is the mismatch still strong (actually stronger!)^. This 

^ We use the same M-R relation in both sets of simulations. 
{Rl + R2)/{R^j^^ + i?!*"') is equal to 4.52 for the realistic stars 
and to 3.30 for the polytropes. Normalising by the half-mass radii 
amounts to a relative contraction of the polytropic models by a 
factor 3.30/4.52 = 0.73. 



fact suggests that it could be possible to scale out much of 
the dependency on the stellar structure by use of some sub- 
tle parameterisation of the "closeness" of the collision that is 
a better representation than dmin/ {Ri + R2) of the amount 
of stellar matter which is highly affected by the collision. In 
cases with stars of very different sizes, a good variable could 
be the mass fraction of the larger star inside dmin or some 
more realistic closest approach distance that includes cor- 
rections for non-Keplerian effects at small distances. In the 
same spirit, rather than using V^'^/V* (or the half-mass ver- 
sion of this quantity) , we could look for a parameterisation of 
the encounter's severity that reflects the energy input com- 
pared to the total binding energy of the stars, for instance. 
In other words, our only hope to find a general description of 
our results that is both relatively simple and robust enough 
to allow some amount of extrapolation, is to trade apparent 
complexity in the results for physically motivated complex- 
ity in the parameters! At any rate clever parameterisations 
can possibly bring together the results of collisions for dif- 
ferent stellar structures only as long as general quantities 
such as the mass and energy losses are concerned. Because 
the entropy and chemical profile of an evolved MS star is 
very different from a homogeneous polytrope, the structure 
and evolution of the collision products strongly depend on 
the use of re^isW^nitia^nod«ls^|^demonsta^e^^^Sil^ 
& Lombardi Tl997l) . 

Such remarks, as well as our comments on the 
strong limitations to the use of published fitting formulae 
(point [(iii)1 above) convinced us that any successful math- 
ematical description of the collisions' outcome should stem 
from physical considerations if it has to be used not only 
as a handy summary of the computed collisions but also to 
extrapolate to somewhat different initial conditions. Unfor- 
tunately, due to the complexity of the physical processes at 
play during collisions, such a "unifying" description seems 
very difficult to find and has escaped us so far. This pushed 
us to cover as completely as possible the relevant domain of 
initial conditions and motivated the use of an interpolation 
algorithm to determine the outcome of any given collision 
with parameters inside this domain. 

3.3 Using the collision results in stellar dynamics 
simulations 

3. 3. 1 The struggle for fitting formulae 

The result of a collision is described through a small set 
of quantities; the fractional mass loss (Mi -I- M2 — Mi' — 
M2')/{M\ -\- M2), the new mass ratio, the fractional loss of 
orbital energy and the angle of deviation of the relative ve- 
locity. Note that these values completely describe the kine- 
matic outcome of a collision only if the centre-of-mass ref- 
erence frame for the resulting star(s) (not including ejected 
gas) is the same as before the collision. Asymmetrical mass 
ejection violates this simplify ing assumption by giving the 
resulting star(s) a global kick iBenz fc HillsllT98^ . However, 
we checked that the kick velocity is generally much lower 
than the relative velocity. Thus, this simplification, which 
greatly reduces the complexity of the situation, should not 
lead to an important bias in the global influence of collisions 
in the energy balance of a star cluster. 

We have kept the final SPH particle configuration for 
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Figure 14. Fractional mass loss as obtained in simulations with polytropes (dashed lines) and realistic stellar models (solid lines). For 
the simulations with polytropes, we used models of indices n = 1.5 and 3 for the small and the large star, respectively. The same M-R 
relation was used for both sets of simulations. In panel (a), we normalised the Keplerian closest approach distance by the sum of the total 
stellar radii. We note that, except for head-on collisions which result in the same mass-loss, polytropes lead to a systematic over-estimate 
of SM. This is probably due of less concentrated density structure of the n = 3 polytrope as compared to a "real" 2 Mq star. In panel 
(b), we use the half-mass radii as a normalisation. In this representation, the agreement is much better up to ^ 2(R^^^ + R2^'')- The 
small bump on the low velocity curve for polytropes at ~ 0.9(-R|^'^' + R'^'') is probably the sign that this collision is a two-stage merger, 
i.e., that a short-lived binary is formed at first periastron passage that merge into a single object at second passage. The symbol type 
indicates the outcome of the collision: triangle for a complete disruption, open square for a binary, filled square for a fly-by and round 
dot when only one star remains (merger or disruption of the smaller star). 



(nearly) all our simulations. This would allow us to re- 
analyse these files and extract other quantities of interest, 
like the amount of rotation imparted by the collision, a quan- 
tity worth investigating because it can deeply influence the 
subsequent evolu tion of the star(s) llMaeder fc Mevnetl200Ct 
ISflls et alJliool and lead to observational signatures that 
would reveal the importan ce of collisions and close en coun- 
ters in given environments (lAlexander fc Kumailliool . An- 
other interesting issue is the resulting internal stellar struc- 
ture. This is key to a prediction of the subsequent evolution 
and obsegatiana^d^tec^^il^^^^olliswi^^ro^^ 
et al. Il997ll200ll for in stance). Unfortunately, according to 
iLombardie^aDlllQQSt) . low resolution and use of particles 
of unequal masses can lead to important spurious particle 
diffusion in SPH simulations so that our models are proba- 
bly not well suited for a study of the amount of coUisional 
mixing, for instance. 

Fig. 1151 shows the (interpolated) fractional mass loss 
in the (dmin, Ke?) plane for various (Mi,M2) pairs. Note 
how the "landscape" changes from one choice of (Mi,M2) 
to another one. Such relatively complex structure obviously 
is a challenge to attempts of describing the results by fitting 
formulae. 

Let's report some unsuccessful attempts at finding easy- 



to-express regularities in the simulation data. We first con- 
vinced ourselves that the outcome of a collision does de- 
pend on both stellar masses and not only on the mass ra- 
tio q = M1/M2. This is demonstrated in Fig llbl in which 
we plot the total fractional mass loss for head-on mergers 
with ~ 0. If this quantity depended on Afi and A'h 
only through q, we would obtain constant S values along 
diagonals, which is not the case. Fig. 1171 depicts another 
wrong guess, namely that for stars of very different sizes, 
the mass loss would only depends on the kinetic energy of 
the impactor (and on dmin) and not on its mass. There is 
not much interest in explaining in detail all the strategies 
we have tried to reduce our huge dataset to a more man- 
ageable mathematical formulation. As a last illustration of 
the difficulty of such a programme, let's mention our at- 
tempt at a global parameterisation of the mass-loss curves. 
We found a 3-parameter formula that allowed good fits of 
individual dmm/(-Ri +-R2) — > 5M/{Mi+M2) relations (for 
fixed Ml, M2 and VJ-S" This looked very promising but we 

^ To achieve these fits, we removed all points corresponding to 
the formation of binaries, because our parameteriscd function was 
monotonically decreasing with increasing dmin E'ud could not re- 
produce extra mass loss due to subsequent periastron passages. 
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Figure 15. Collisional fractional mass losses for four different {Mi, M2) pairs (values in Mq). Simulation Data. White dots show the SPH 
simulations. The contours and colour maps are a bi-cubic interpolation of the SPH results. S is the fractional mass loss S = SM/ (Mi + M2). 
Masses are in units of Mq. In each frame, the upper left contour indicates fractional mass loss larger than 85%. 



were left with 1180 sets of parameters to be adjusted in turn 
by some "meta-formula" , with the added difficulty that they 
displayed a lesser level of regularity than the raw collisional 
data itself. This proved unmanageable. Furthermore, this 
parameterisation had no sound physical justification. 

To end this subsection on a more positive note, let's 
turn to Fig. 1181 In this diagram, we plot the relative mass 
gain or loss for the larger star, S2, as a, function of the usual 
half-mass normalised A and v. The figures are remarkably 
smooth in the sense that collisions with comparable mass 
ratios, but otherwise different Mi and M2, and same (A, v), 
produce very similar 32- There is thus some hope that, in 
further investigations, we could discover some "universal" 



^2 = 52{q,i',\) relation to describe this regularity. Such a 
description would be particularly useful to explore, with ana- 
lytical or semi-analytical models, the possibility of run-away 
merging sequences in the evolutions of dense clusters. Using 
realistic SPH results to re-examine these scenarios is one 
important application of the present work (see Sec. II. It . 



3.3.2 Interpolation of the collision results 

Being unable to distillate the results of our SPH simula- 
tions into any compact mathematical formulation without 
losing most of the information, we resorted to the following 
interpolation strategy. In the 4D initial parameter space. 
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Figure 18. Plots of the relative modification of the mass M2 of the larger colliding star as a function of dn 
colour-coded in red to yellow, are normalised as fractions of A/2. Mass increases, colour-coded in green tones are normalised as fractions of 
Ml. We choose this two different normalisations so these relative mass changes are always comprised between and 1 in absolute value. 
(Mj — M2)/M2 can me interpreted as the "fractional damage" caused by the "bullet" (small star) to the target, while (Mj — M2)/Mi 
is the "efficiency" by which mass of the small star is added to the more massive one. 



6(d . =0, V ,) l.OQ^V JV^^1.Q4 

^ min ' cont/ cont' * 



10= 



10 



0.1 



mil 
I IMF 



1 3.7x 
1 9.7X 
1 2.3x 
I 5.4X 
1 l.Ox 
1 1.5x 
1 2.1X 
1 3-Ox 
I 4. Ox 
I 4 

I 5.4X 
I 6.2x 
1 7. Ox 
1 7.4x 
1 7.9x 



10-*<a£9 
10-*<SS2 
10-3<iS5 
10-=<S£1 
10-=<S£1 
10-2<(5S2 
10-2<5S3 
ia-^<iS4 
10-=<as4 
10-'<6S5 
10-2<S£6 

io-2<as7 

10-2<iiS7 
10-=<5fi7 

io-2<s£a 



0.1 



10 



Mj in 



7x10-* 
3xl0-» 
4x10-1 
0x10-2 
5x10-2 
1x10-2 
0x10-2 
0x10-2 
6x10-2 
4x10-2 
2x10-2 
0x10-2 
4x10-2 
9x10-2 
3x10-2 



10^ 



Figure 16. Fractional mass loss (5) for all head-on, "zero- 
velocity" collisions. Here, Vcont is the contact velocity (at sep- 
aration Ri + R2)- S is not constant on lines of constant M1/M2. 



the simulations form a irregular grid of points. We com- 
pute a Delaunay triangulation of this set iSedgewick|[i988l 
Chap. 28) using the program QHULL* iBarber et alJll99el) 
which allows us to interpolate the results onto a regular 4D 
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Figure 17. Fractional mass loss for all collisions between a star 
of mass 19.3 M0 and stars that are at least 10 times less massive. 
Here we normalise the distance by R2 , the total radius of the large 
star. Points are colour-coded according to the kinetic energy of 
the impactor at contact (in solar units, GMq/Rq). 



grid. To evaluate the value of any of the four quantities that 
summarise the outcome of a collision, £5, we first find the 
simplex S of the triangulation, if any, that contains the 4D 
point P of the initial conditions of the collision. This sim- 
plex has 5 vertices: Qi,i = 1 . . .5. By removing one of these 
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Figure 20. Four slices through our interpolation grid for the coUisional fractional mass loss. Wc performed cuts that correspond to 
(Mi,M2) values close to those of Fig. 1151 <5 is the fractional mass loss 5 = <5M/(Mi + Ml). Masses are in units of Mq. The Delaunay 
interpolation produces some artifacts at low and high relative velocities, in particular for the top-left and bottom-right panels (compare 
with Fig. 1151 . This is due to the simplices being very elongated near the border of the convex hull of our data (in four dimensional space). 



vertices, say Qk and replacing it by P, one forms another, 
smaller simplex Sh that is contained in S. We compute the 
interpolated value of Q at point P, Q(P) from its values at 
the vertices Qi, 0(Qi) by linear combination with weights 
yst^l^s where Vs stands for the (hyper-)volume of <S. Of 
course, this procedure does not allow extrapolation outside 
the convex hull of our simulation initial conditions. Another, 
more tricky problem is that, near the borders of this convex 
hull, simplices can be very elongated which means that the 
interpolations can be done with data points corresponding 
to very remote initial conditions instead of using more lo- 



cal information. This is illustrated, for two dimensions, in 
Fig. EH 

However, we did not find a better procedure. We tried 
to use a kernel-based method, a la SPH, but it produced 
very poor results. The main problem with this class of al- 
gorithms is to adapt locally the 4 independent axis of the 
kernel ellipsoid in such a way that only neighbouring data 
points contribute to the evaluation at a given point. Defining 
"neighbours" is unfortunately not obvious in a parameter 
space with no natural metrics^. 



^ This problem about the metrics being ill-defined actually 
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Figure 19. Dclaunay triangulation in two dimensions. In our 
interpolation method we would get the value for the point P 
(small square) from data points 1, 2 and 3 that are the vertices of 
the triangle P lays in. Although they are much closer to P than 
1 or 3, points 4 and 5 would not contribute at all! 



The quality of the data obtained through our interpo- 
lation mechanism is illustrated on Fig. 1201 It shows four 2D 
slices of the fractional mass loss interpolated onto the 4D 
grid. Each slice corresponds to a (Mi, M2) pair chosen as to 
be as close as possible to the values used for Fig. 1151 allow- 
ing a direct comparison. The general dependency of the mass 
loss on impact parameter and relative velocity is well repro- 
duced but some details are smoothed out while small arti- 
facts have appeared near the borders of the domains for the 
reason explained above. We interpret the horizontal "penin- 
sula" of high mass loss for Mi = M2 = 0.5 Mq as the result 
of the interpolation between a simulation which was inte- 
grated long enough for complete merger (visible on Fig.lTsl 
for the lowest relative velocity value at dmin/ {R1+R2) = 0.8) 
and led to relatively high mass loss and others which were 
stopped at an earlier phase (unmerged binary). 

The table thus computed is the backbone of the rou- 
tine that implements stellar collisions in our Monte Carlo 
simulations of stellar clusters. Collision outcome quantities 
are indeed easily obtained through a second, much quicker, 
interpolation stage using this regular grid. Of course, extrap- 
olation prescriptions have to be specified for events whose 
initial conditions fall outside the convex hull of the SPH sim- 
ulation points. Most commonly, this happens when a col- 
lisionally produced star with mass outside the 0.1-74 M0 
range experience a further collision. In such cases, we try to 
re-scale both masses while preserving Mi /M2 to get a "sur- 
rogate collision" lying in the domain covered by the SPH 
simulations. If Kci is too low or too high, we increase or de- 
crease it to enter the simulation domain^". In its present 
state, this treatment of "extreme" velocities is not com- 



also plagues the Dclaunay method, but in a less visible and 
apparently less harmful way! To get more acquainted with 
Dclaunay triangulation, see the interactive demonstration at 
[http : //wwwpl6 . f ernuni-hagen.de/GeomLab/VoroGlide/ 

All this fiddling does not violate mass or energy conservation 
as collision results are coded in a dimensionless fashion in the 
interpolation grid and are scaled back to the real physical masses 
and velocities before they are applied to the particles in a stellar 
dynamics simulation. 



pletely satisfying. At very high Kei, our data do not show 
convergence toward a unique mass loss curve. Instead, the 
domain of complete disruption keeps extending to higher 
and higher impact parameters with a progressive steepen- 
ing of the mass loss curves for "fly-bys" . At very small ve- 
locities, the values of the table can be trusted only when 
no binary has formed or if the binary evolution has been 
followed up to merging. In case of binary formation, some 
constant fractional mass loss could be used in order to re- 
flect our finding that the process of binary merging, that 
requires more and more pericentre passages for larger and 
larger dmin, eventually leads to an amount of mass loss rela- 
tively independent of this impact parameter. A more precise 
determination of the maximal dmin that still leads to binary 
formation for small initial velocities would allow us to know 
where to switch between this prescription and interpolation 
in the table. Finally, cases with too high dmin are treated as 
purely Keplerian hyperbolic encounters with no mass loss, 
which is a very good approximation. Recently, in the frame 
of our work on coUisional run-away formation of very mas- 
sive stars, we have implemented a few more small tricks to 
complem ent our "blind" interpo lation routine and reduce its 
artifacts llFreitag et al.ll2005bla|) . 

An important aspect of the work reported here is that 
we make the data describing the initial conditions and 
outcome of all our simulations available on the web, on the 
site of the "MODEST" working group on stellar collisions at 
|http : //obsww . unige . ch/~f reita g/MDDEST_WG4/FB_Col lision_Data/| 
We provide a description of the outcome of a given collision 
in terms of the number and masses of star(s) at the end of 
the simulation and their orbital properties. Colleagues are 
invited to develop their own methods to use this data and 
share their experience with others, including the authors of 
the present paper. Files containing detailed information for 
all SPH particles at the end of a simulation are available 
upon request to MF. 



4 CONCLUSIONS AND FUTURE WORK 

In this article, we presented a large set of simulations of col- 
lisions between two main sequence stars. More than 14 000 
SPH simulations have been computed over about four years 
to complete that database. Initial conditions span M, — 
0.1 — 75 M0 for the stellar masses, impact parameters cor- 
responding approximately to dinin/(i?i + R2) = — 0.9 and 
relative velocities at infinity ranging, more or less, from 0.03 
to 30 times the stellar escape velocity. This represents a ef- 
fort of unprecedented breadth in this field. 

Our motivation in this work was to incorporate the ef- 
fects of stellar collisions into models of dense stellar systems 
like galactic nuclei with as much realism as possible. To reach 
this goal we developed a module that interpolates between 
our results to predict the outcome of any collision with ini- 
tial conditions inside the (large) domain of parameter space 
we explored. Results of our dynamical simulations of dense 
clusters including collisions are Dresented elsewhere (Freitag 
fc Benz 2001a., 2002.; Rasio ct al.,,2004: .Freitag ot al...2004bja[ 
l2005bla[K 

The quest for a handy mathematical description of these 
results has been unsuccessful so far. This was a source of 
disappointment but we hope that further study of our sim- 
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ulation data will eventually reveal some way of casting our 
results in a compact and physically enlightening formula- 
tion. Exploring when and why the excellent match with the 
SS66 prescription is found is a possible way to this goal. 

Beyond the scientific exploitation of this important 
database, -either to develop a better understanding of the 
physics of collisions or as an ingredient for coUisional stellar 
dynamics studies- we are aware that other kinds of stellar 
collisions, not treated here, are also of great astrophysical 
interest. First, in galactic nuclei, collisions with red giants 
are certainly more frequent than MS-MS encounters. So we 
should assess the importance of this process (for stellar evo- 
lution and stellar dynamics) and, if needed, compute a set 
of simulations similar to the one presented here. This would 
complement the work of iBailev fc Davie3 ^1999^ . Collisions 
with compact remnants should also be taken into account 
as they may be of great importance in producing peculiar 
objects. 

Our work is not well adapted to the physical conditions 
that are typical in globular clusters. Indeed, we did not study 
with particular care the low velocity, quasi parabolic encoun- 
ters, by far the most common in those environments. How- 
ever, they have already been quite thoroughly studied by 
other researchers. The conditions for tidal binary formation 
are probably better determined using other methods and 
their lon g term evolutiori, w hose nature and result is still 
debated iMardlinelll995albl) . extends over much too many 
hydrodynamical time scales to be tackled by SPH. Still in 
the domain of globular clusters, the study of hydrodynam- 
ical efl^ects (including direct collisions) in interactions with 
or between binarv stars is still in its infanc^^Goodman & 
Hernquist 1991; Davics ct al 1993. .1994) . 

Finally, going back to galactic nuclei in which we are 
mostly interested, let's mention that the connected problem 
of the destructive close encounter between a central massive 
black hole and a star, even though it has been the focus of 
many papers JCarter fc Luminedll983t iBicknell fc Gingoldl 



198alEvans fc Kochanekll989HLaguna et alll993HFulbrigh1 



199( ; iMarck et aLlll996t lAval et alJl2000l: iBogdanovic et al 



20041 and many others), has still not been explored sys- 
tematically. For instance, all studies so far have assumed 
simplified stellar models. In this problem, however the main 
uncertainties probably lurk in the post-disruption evolution 
of the stellar gas, rather than in the "collision" process itself. 



APPENDIX A: SMOOTHED PARTICLE 
HYDRODYNAMICS STELLAR MODELS. 

Al Stellar structure used in the simulations 

In Fig. lAll we show the density profiles for some realis- 
tic models for MS stars and compare them with polytropic 
stars. It is readily seen that for M > 0.4 M©, stellar struc- 
tures are highly non-homologous and that polytropes do not 
match, even if allowance is made for a Af-variable n index. 
If such a fit is required anyway, a value of n ~ 3.5 seems 
more appropriate for M ^ 1 Mq than the commonly used 
n = 3. 

An SPH particle configuration for a star is illustrated 
on Fig. IA2I One sees that the outermost layers of the star 
are very poorly modelled, with a clear failure at precisely 



reproducing the real stellar radius. Fig. IA3I is a compari- 
son between the density and internal energy profiles of two 
stellar models and their SPH approximation for increasing 
number of particles. 



A2 Choice of the particle number 

In Fig. m we show how the overall results of SPH colli- 
sion simulations (mass and energy loss, deflection angle) de- 
pend on the resolution, i.e. the number of particles used to 
represent the stars. We considered resolution ranging from 
1000-f 2000 to 2000+32 000 particles. 



APPENDIX B: RESULTS OF SPH 
SIMULATIONS 

Bl A few specific colhsion simulations 

Precise descriptions of the physical mechanisms at play dur- 
ing stellar colhgion^ia^^^rea^vbegi^21lMiSij£^^S£ii2_,^ 
Hills .1987. .19921 iLai et alJll993t iLombardi et alJll99d) . 

In this subsection, we just highlight a few particular col- 
lisions from our survey for illustration purposes. We do not 
particularly concentrate on "classical" typical cases because 
they have been well covered in these previous works. In- 
stead, we concentrate on simulations with parameters lying 
on the border-lines of the various regimes. Many of them 
have been re-computed in order to test surprising results. 
Indeed, for lack of sufficient disk space for data storage, only 
the final "state" of each SPH simulation was conserved for 
most runs. So, when any doubtful result appeared, we had 
to re-compute the complete simulation and write the data 
to disk frequently in order to understand the evolution of 
the system. 

In Fig. IBll we show an ofi'-axis low velocity encounter 
between identical stars. As the impact parameter is small, 
the stars merge together after the first periastron passage. 
The colour mapping used in these diagrams allow to trace 
each particle back to its initial radial position in the colliding 
stars. Despite the rather low resolution (about 8000 parti- 
cles in total) a tighter and tighter spiral pattern is clearly 
visible. As explained by Lombardi and collaborators, in low 
velocity collisions, specific entropy s is nearly conserved as 
shock heating is weak, and, as stability of the resulting star 
imposes ds/dr > 0, low entropy material that was initially 
at the centre of the stars, settles in the centre of the merged 
object. In fact, for such gentle encounters, one can even 
predicts the final material stratification by sorting mass el- 
ements from the two stars i n ord er of increasing entropy 
jLombardi et alJll996l l200l l2003h . A consequence of this 
mechanism is that, in mergers between stars with unequal 
masses, the core of the smaller one, having the lowest en- 
tropy, sinks to the centre of the merger. This is also what 
happens in the two collisions depicted in Figures IB2I IB3I 
and IB4I The second collision is an example of a high ve- 
locity merger which produce an object with a total mass 
slightly lower than those of the initial larger star. Such a 
case lies in the tip of merger region in a diagram like those 
in Fig. 9 of the main paper. 

In Figures IB5I and IB6I we display snapshots from one 
of the few head-on collisions in which the small star pass 
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Figure Al. Density profiles for realistic star models iSchaller et al.lll992t ICharbonnel et aljITooih for low (top) and high (bottom) 
mass stars. The dashed lines are polytropic models for n = 1.5, 3 and 4, in order of increasing concentration. Below 0.4 Mq, the density 
structure is well represented by a polytrope with n = 1.5 but no good polytropic fit is possible for higher masses. 
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Figure A2. SPH realisation of a SMq 
stellar model with ~2000 particles. Round 
dots show the positions of SPH parti- 
cles with a symbol surface proportional to 
the particle's mass. The big dashed cir- 
cle shows the size of the star according to 
the structure model. Plain line circles de- 
pict the concentric spheres on the surface 
of which the particles' centres are placed. 
Particles on the > corner of the di- 
agram have been removed to show the ac- 
tual half-size of particles on each sphere 
(dotted circles of radius h). 
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Figure A3. Comparison between the theoretical internal structure of a 2 Mq stellar model (solid lines) and SPH realisations of it with 
increasing number of particles (dashed and dotted lines). We show the mass density (left column) and the internal energy (right column). 
In the top row, we use the radius as the abscissa, while we use the enclosed mass for the diagrams of the bottom row. The top plots 
show that the outermost part of the star are poorly represented. However, it is clear from the bottom plots that only a tiny mass fraction 
suffers from this mismatch. The SPH profiles are the kernel-interpolated values along the line z = 0. 



through the large one and remains essentially intact. A fur- 
ther peculiarity of this colhsion is its relatively large mass 
ratio: q = 0.24. No collision with a larger q resulted in a 
similar outcome. Figures inn and depict a more typical 
"fly-by" in the sense that it has non-vanishing impact pa- 
rameter. However, this high velocity encounter lies very close 
to the strip of complete disruption of the small star. For this 
particular simulation, the small star loses more than 89 % of 
its mass! The remaining cloud has a very low central density, 
around 10~*gcm~^. It is made of only 187 particles so sim- 



ulations with higher resolution are clearly needed to confirm 
that the production of such tiny survivors is not a numerical 
artifact. It is however unlikely that such small, rarely formed 
objects, may have important astrophysical relevance, either 
as detectable "exotic" stars or dynamically. 

We finally present a collision from the high velocity, 1- 
star branch, i.e. a case of collisionally induced evaporation 
of the small star. We particularly checked that this kind of 
outcome was real and not some artifact cause by our anal- 
ysis software. Indeed, during these verifications, we noted 
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Figure A4. Study of the dependency of collision results on the particle number for 4 sets of collision simulations, (a) Fractional mass 
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(<5i?orlj /^orb = !)■ (c) Deviation of the coUisional deflection angle ficoll from the value for pure Keplerian point-mass trajectory Sgrav 



that many cases of nearly complete destruction of the small 
star Uke the one described in Figures IB5I and IB6I were mis- 
interpreted because our code missed the second, much lower, 
density peak. Consequently, we had to re-analyse all high v 
collisions that were reported to result in the disruption of 
the small star. We conclude that although the precise loca- 
tion of its right (large dmin) edge may depend on numerical 
issues (resolution, analysis procedure), the 1-star branch is 
real. Inspecting the last frames of Figures IB9I and IB 1 01 and 
Fig. lBlll makes this fact obvio us. Furth ermore, such collision 
results have been reported by iLai et al.i (^IflflAl . 

B2 Examples of mass and energy loss results 

Fig. IB12I shows the energy and mass loss curves for the 
simulations of collisions between stars with (Mi, M2) = 
(0.5, 12)M0, and (12, 12)Mo. Similar curves for other 
(Ml, M2) couples are available upon request to MF. The 
can also easily be drawn using the complete tables of colh- 
sion results available on-line. 
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Figure Bl. Collision between two 2.0 stars at V^ylV* = 0.77 (465 kms^^) and dmin/C^?! +^^2) = 0.1. Each plot show the position, 
projected onto the orbital plane, of SPH particles that lie close to this plane. Beware that the length scale may change from frame to 
frame. This collision creates a merged star with mass 3.81 Mq. The particles arc colour-coded according to their rank in the initial stellar 
models. 
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Figure B2. Collision between stars of masses 1.0 and 3.0 at V^^/K = 0.07 {44kms~i) and dmin/C-Ri + R2) = 0.0. This collision 
creates a merged star with mass 3.80 Mq. The particles are colour-coded according to their rank in the initial stellar models. 
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Figure B3. Collision between stars of masses OAMq and 4.0 Mq at V^'^JVt = 3.24 (2180 kms"!) and (imin/{^?i -f = 0.04. The 
colours code the gas density. This collision results in a merged star with mass 3.75 Mq. 
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Figure B4. Snapshot of the collision simulation of Fig. IB3l at later stage. One sees the small star spiralling into the centre of the larger 
star. Panel (a); density plot. Panel (b): Velocity plot with colours indicating the radial position of each particle in the initial stellar 
models. Particles from the large star are coded in red to yellow, from centre to surface. Particles from the small star are coded in dark 
to light green. The small dots are the position of the particles. The velocity scale is given by the horizontal line segment at the bottom 
right of the diagram. 
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Figure B5. Collision between stars of masses 1.7 Mq and 7.0Mq at V^^JV, = 3.68 (2620 kms"!) and dmin/{.R.i -l-^?.2) = 0.0. Not only 
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Figure B7. Collision between stars of masses 0.5 M© and 2.OM0 at V^^^/V* = 4.48 (2620 kms^^) and dmin/C-Ri + R2) = 0.15. 
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Figure B9. Collision between stars of masses 0.4 A/q and 1.7Mq at V^^/V* = 3.80 (2180 kms"!) and dmin/(^f.l + -^2) = 0.11. 



40 M. Freitag and W. Benz 



T= 1.50x10' hours 




p in M^Rg' 

55.36x10-'' 
59.75x10-'' 
5 1.73x10-* 

Ja.esxio" 

S 4.4ixl0-' 

Se.esstio-" 

J 1.04j<1G-s 
5 2.33x10-* 
5 7.04x10-* 
5 9.07x10-* 
2 8.55x10-" 
56-91x10-= 
53.60x10-" 

~ i.iaxio« 



T= 1.11x102 hours 



5 
X in Rn [xlQi] 



o 

X 



a; 



-4 - 



T 



i 



' _- A. 



_L 



_L 



X in 



2 4 



P in MgRg= 

— S.OSJtlQ-" 

Je.-isxio-* 
5e.e6xio-» 

5 1.35x10-" 
5 £.15x10-" 
S 3.70x10 " 
! 8.54X10 » 

Sa.saxio-'^ 
5 i.aixio-' 

5 3.33x1 Q-s 

5e.o5xiQ-' 

55.74x10-5 

Sa.osxio-' 

5 1.36x10-1 
54.30x10-1 
5 1.08x10* 
" 3.06x10" 



Figure BIO. Continuation of the sequence of Fig. IB9 
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Figure Bll. Enlargements of the last frame of Fig. IBlOl Position and velocities are relative to the particle with the highest density 
in each panel. Panel (a): surviving core of the large star, a 1.26 Mq rotating star. Panel (b): remaining of the core of the small star, 
an unbound expanding gas cloud. The velocity of this cloud in the centre-of-mass reference frame of the collision is nearly 1000 ism s^^ 
while the velocity of the small star was initially 17701tms~^. 
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Figure B12. Relative mass and kinetic energy losses for all collision simulations between stars of masses 0.5 and 12 Mq (column (a)) 
and 12 and 12 (column (b)). Half-mass radii arc used to normalise parameters. Tcont is the orbital kinetic energy at "half-mass 



contact" (separation equal to -\-R2^'), assuming purely Keplerian acceleration. For very small Vj.^, all encounters result either in 
mergers or in bound binaries that should eventually merge together, with fSTcont /Tcont = 1 and a higher SM/ (Mi -|- M2) as consequences. 
At high velocities, the domain of 100 % complete energy loss extends to dj^iu values where mass loss is only partial. This is due to the 
complete disruption of the smaller star after it emerges from the large one. The shoulders on the low velocity mass loss curves are due 
to the formation and subsequent merging of a binary. 
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